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1.0  Introduction 


This  document  comprises  Part  1  of  the  final  report  for  the  research  effort  entitled  “GPR 
performance  in  the  presence  of  buried  biomass”  (Contract  DAAD 19-01- 1-0760,  Project  No.:  E- 
43042-EL-000-0 1234-1).  There  is  extensive  anecdotal  evidence  that  biomass  (i.e.,  plant  and  tree 
roots)  adversely  affects  ground-penetrating  radar  (GPR)  performance.  In  this  effort  the  Ohio 
State  University  (OSU)  ElectroScience  Laboratory  (ESL)  described  the  physics  that  govern  the 
propagation  and  scattering  of  GPR  waveforms  in  the  presence  of  buried  biomass.  This  work  is 
being  performed  in  collaboration  with  the  US  Army  Research  Laboratory. 

This  report  is  presented  in  two  parts.  In  Part  1  we  document  the  results  of  experimental  GPR 
measurements  performed  at  Site  C62  at  Eglin,  AFB,  FL.  In  Part  2  we  describe  numerical 
modeling  work  that  permits  us  to  simulate  and  interpret  many  of  the  results  presented  herein. 

This  document  is  organized  in  three  major  sections.  In  Section  2  we  discuss  measurements 
performed  with  a  ground-contacting  GPR  antenna.  In  Section  3  we  discuss  surface-scattering 
measurements  performed  with  an  elevated  antenna.  Finally,  in  Section  4  we  describe 
measurements  of  some  supplemental  data,  including  descriptions  of  the  root  structures  found  at 
the  test  site.  Concluding  remarks  appear  in  Section  5. 
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2.0  GPR  Measurements  with  a  Ground-Contacting  Antenna 


A  team  from  ESL  visited  Site  C62  at  Eglin  AFB,  FL  during  December  10-12,  2001  to  perform 
radar  measurements  and  to  collect  samples  of  soil  and  roots.  The  area  is  covered  with  tufts  of 
grass  and  short  scrub  oak  trees  that  are  periodically  cut  back  to  ground  level  (see  Figure  1).  At 
the  time  of  these  measurements,  the  average  height  of  those  trees  was  less  than  24  inches.  Using 
the  local  density  of  trees  as  a  guide,  it  was  postulated  that  the  regions  in  Figure  1(a)  and  (b) 
would  have  relatively  high  and  low  root  densities  respectively.  All  of  the  results  documented  in 
this  report  were  acquired  at  these  locations,  which  are  separated  from  one  another  by  only  a  few 
tens  of  meters  at  Site  C62.  At  each  location  we  selected  regions  of  size  two  meters  square  for  our 
study.  The  region  with  high  root  density  is  referred  to  hereafter  as  “Site  1”  and  the  low  root 
density  region  is  designated  “Site  2”. 


(a)  Site  1  (high  root  density) 


(b)  Site  2  (low  root  density) 


Figure  1.  Test  locations  at  Eglin  AFB,  Site  C62 


2. 1  GPR  Description 

An  ultra-wide  bandwidth,  multi-polarization,  ground  penetrating  radar  (GPR)  was  used  to  collect 
data  at  Sites  1  and  2.  The  “radar”  comprises  a  commercial  RF  vector  network  analyzer 
(HP8753C)  for  collecting  reflection  data  (via  the  S  parameters  Snand  S22)  and  transmission  data 
(S21)  and  a  dual-polarized  antenna  utilizing  a  dielectric-loaded  horn  design  developed  at 
OSU/ESF.  The  horn  section  is  filled  with  an  artificial  dielectric  with  a  relative  dielectric 
constant  of  five,  which  is  well  matched  to  the  local  soil.  The  complete  radar  and  a  conceptual 
view  of  the  antenna  are  shown  in  Figure  2.  The  antenna  was  oriented  such  that  one  of  the  two 
dipole-elements  was  parallel  to  the  direction  of  scan.  The  return  from  the  parallel  arm  comprises 
channel  2  and  the  return  from  the  perpendicular  arm  comprises  channel  1.  Data  were  collected 
over  the  frequency  range  0.05  to  2.0  GHz. 
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Conductor 


(a)  horn  antenna  and  cart-mounted  GPR  (b)  horn  antenna  schematic 

Figure  2.  Ground-contacting  GPR  system  with  dielectric-loaded  horn  antenna.  The  horn  is  made  from  two  folded 

bow-tie  elements  terminated  in  resistive  film. 


The  radar’s  control  and  data  collection  functions  were  performed  by  a  laptop  computer  using 
OSU-developed  software  operating  over  a  GPIB  connection.  The  configuration  of  the 
components  is  illustrated  in  Figure  3  and  Figure  4.  The  network  analyzer  and  laptop  were 
mounted  on  a  non-metallic  wheeled  platform.  The  antenna  was  placed  on  a  circular  sled 
attached  to  the  back  of  the  wagon  by  two  trailing  arms. 

Data  were  acquired  using  this  system  at  41  positions  spaced  at  three-inch  intervals.  Herein  we 
refer  to  this  spatial  sampling  as  a  “scan”.  At  each  position  in  the  scan  three  polarimetric 
components  were  sampled  as  Sn,  S22  and  S21.  Each  sample  comprised  measurements  at  401 
frequency  steps  from  0.05  to  2.0  GHz  (step  size=4.875  MHz).  The  IF  bandwidth  was  set  to  300 
Hz,  and  the  power  level  was  set  to  0  dBm.  The  real  and  imaginary  parts  of  the  measured  S 
parameters  were  stored  as  binary  files  using  the  format  described  in  the  Appendix. 


Figure  3.  GPR  system  (top  view).  A  small  wheeled  platform  is  used  to  drag  the  ground-contacting  antenna,  which 

is  mounted  on  a  circular  plastic  sled. 
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Figure  4.  GPR  system  (side  view).  The  antenna  is  connected  to  the  network  analyzer  via  a  semi-rigid  coaxial 
cable.  Trailing  arms  connect  the  antenna  platform  and  the  wheeled  cart. 


Prior  to  making  any  measurements,  we  removed  the  large  surface  vegetation  (primarily 
the  low  scrub  oak  trees).  Local  EOD  personnel  then  performed  a  sweep  of  the  areas  using 
magnetometers  and  metal  detectors.  (Site  C62  is  known  to  contain  unexploded  ordnance.)  At 
Site  1  no  foreign  objects  were  found,  but  at  Site  2  a  length  of  wire  and  some  metal  fragments 
were  located  and  removed.  Site  2  also  contained  a  segment  of  steel  rebar  that  had  been  buried  in 
a  vertical  orientation  for  use  as  a  fiducial  marker.  That  rebar  was  temporarily  removed  during 
our  measurements  and  later  returned  to  its  original  location. 

We  next  dug  an  access  trench  to  permit  a  steel  rod  to  be  inserted  into  each  test  area  as 
shown  in  Figure  5.  The  target  rod  was  3  meters  long  with  a  0.875-inch  diameter.  The  access 
trench  was  located  adjacent  to  the  midline  of  each  test  area  and  had  a  depth  of  approximately  5 
feet. 

For  each  of  the  two  sites,  seven  separate  scans  were  performed  as  indicated  in  Table  1. 
The  first  scan  was  performed  prior  to  inserting  the  test  rod.  The  radar  path  was  perpendicular  to 
the  axis  of  the  (absent)  rod.  After  this  initial  work,  the  steel  rod  was  inserted  at  a  depth  of  three 
feet.  With  the  rod  inserted,  the  GPR  was  scanned  over  the  rod  twice:  first  through  the  centerline 
of  the  rod  and  a  second  time  offset  from  the  center  away  from  the  access  trench.  Hereafter  we 
refer  to  these  scans  as  “centered”  and  “offset”,  respectively.  The  offset  path  is  offset 
approximately  three  feet  from  the  rod’s  centerline  in  the  direction  away  from  the  trench.  The  rod 
was  then  removed  and  re-inserted  at  a  depth  of  two  feet.  The  centered  and  offset  scans  were 
repeated.  Finally,  the  centered  and  offset  scans  were  repeated  for  the  rod  at  one  foot  depth. 
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(Note  the  use  of  decreasing  target  depth  to  avoid  disturbing  the  soil  between  the  radar  and  the 
target.) 


Table  1  GPR  measurement  scans  taken  at  Eglin  AFB. 


Center  Measurement  Pass 

Offset  Measurement  Pass 

Rod  at  one-foot  depth 

Rod  at  one-foot  depth 

Rod  at  two-foot  depth 

Rod  at  two-foot  depth 

Rod  at  three-foot  depth 

Rod  at  three-foot  depth 

No  Rod 

(a)  Side  view  (b)  Plan  view 

Figure  5  GPR  measurement  paths  (not  to  scale). 


2.2  Results 

The  GPR  data  was  processed  using  software  developed  in  Matlab  to  remove  the  “background” 
response  (i.e.,  reflections  from  the  antenna,  cabling,  and  other  components  in  the  system).  This 
was  accomplished  in  the  frequency  domain  by  computing  the  average  response  over  all  41 
spatial  positions,  and  subtracting  that  response  from  the  entire  scan.  Examples  of  this 
“background- subtracted”  data  are  shown  in  Figure  6,  Figure  7,  and  Figure  8  for  a  scan  over  the 
center  of  the  one-foot-depth  rod  at  Site  2.  A  large  amount  of  high  frequency  (>1.5  GHz)  clutter 
appears  in  the  S22  channel  (see  Figure  7)  and,  as  noted  above,  the  frequency  domain  data  were 
truncated  at  1.5  GHz  to  eliminate  this  clutter. 
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Figure  8  S2i  Site  2,  frequency-domain  one-foot-depth  rod,  centered 


The  complete  set  of  time-domain  responses  for  Sites  1  and  2  are  shown  in  Figure  9  through 
Figure  50  for  the  frequency  band  50-1500  MHz.  The  response  of  the  rod  at  three  feet  is  shown  in 
Figure  9  through  Figure  20,  the  two-foot  rod  results  appear  in  Figure  21  through  Figure  32,  and 
the  one-foot  rod  results  are  found  in  Figure  33  through  Figure  44.  The  response  of  the  empty 
sites  are  presented  in  Figure  45  through  Figure  50.  The  same  amplitude  scale  is  used  in  all 
figures. 

The  time  of  the  surface  return  is  approximately  -6.5  ns  in  these  Figures.  To  suppress  strong 
near-surface  returns,  the  first  6.5  ns  of  the  time  domain  returns  have  been  attenuated  by  a  factor 
that  varies  linearly  from  -40  dB  to  0  dB  over  the  time  interval  0  to  -6.5  ns.  (Data  after  6.5  ns  are 
not  affected  by  this  processing.)  The  soil  dielectric  constant  at  this  location  is  roughly  five  (see 
Part  2  of  this  report),  which  implies  a  round-trip  propagation  velocity  of  2.6  inches  per 
nanosecond.  From  subsequent  excavations  of  the  site  (described  in  Section  4),  we  know  that  the 
roots  reside  primarily  at  depths  of  less  than  12  inches,  which  correspond  to  times  of  -6.5  to  -11.1 
ns. 
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The  position  of  the  rod  is  indicated  by  arrows  in  Figure  9  through  Figure  12.  The  Site  1,  center- 
scan  data  contain  an  artifact  (a  strong  ringing  that  persists  at  all  depths)  near  the  -35  inch 
position.  The  source  of  that  artifact  is  unknown,  but  is  probably  the  result  of  scattering  from  a 
small  piece  of  metal  (e.g.,  a  dropped  screw)  buried  just  below  the  surface  at  that  location. 

As  an  example  of  these  data,  consider  the  response  of  the  rod  at  three-foot  depths.  The  relevant 
data  for  Site  1  are  found  in  Figure  9,  Figure  13,  and  Figure  17.  The  arcs  that  appear  in  the  time 
interval  -20  to  -30  ns  at  positions  -30  to  +40  inches  are  the  back- scattered  target  response,  which 
have  the  hyperbolic  shape  characteristic  of  a  point  scatterer.  The  S22  channel  has  the  strongest 
response,  since  that  antenna  arm  is  co-polarized  with  the  target  rod.  The  response  of  the  cross- 
polarized  (S21)  channel  is  relatively  weak.  As  expected,  rods  at  shallower  depths  produce 
stronger  responses. 

The  response  of  the  roots  appears  to  be  weak  and  not  easily  characterized.  There  is  little 
difference  between  the  early-time  (-6  to  -10.6  ns)  responses  at  Sites  1  and  2.  The  Site  1  data  do 
show  somewhat  more  scattering  in  the  region  occupied  by  the  roots,  but  there  are  no  clear 
isolated  point-like  scatterers.  The  early-time  response  at  Site  1  is  comparable  for  both 
polarizations,  as  one  would  expect  for  randomly  oriented  targets.  Since  early-time  responses  are 
also  seen  at  Site  2,  those  returns  might  be  attributed  to  antenna  mismatch  caused  by  the  rough 
ground  (i.e.,  randomly  placed  tufts  of  grass)  or  inhomgeneities  in  soil.  The  modeling  analysis 
presented  in  Part  2  will  offer  evidence  that  these  returns  are  due  to  the  roots. 

The  cross-polarized  (S21)  return  offers  perhaps  the  best  visible  indication  of  roots.  The  results 
show  that  while  this  component  is  weak,  the  cross-polarized  early-time  response  for  Site  2  is 
consistently  weaker  than  the  response  from  Site  1.  This  finding  has  a  plausible  physical 
explanation:  since  roots  have  a  specific  but  random  orientation,  roots  would  tend  to  produce  a 
random,  but  non- zero  cross-polarized  response. 


Figure  9  Su  Site  1,  center  scan  of  three-foot-depth 
rod  . 


Figure  10  Su  Site  1,  offset  scan  of  three-foot-depth 
rod  . 
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Figure  11  Su  Site  2,  center  scan  of  three-foot-depth 
rod. 


Figure  12  Su  Site  2,  offset  scan  of  three-foot-depth 
rod. 


Figure  13  S22  Site  1,  center  scan  of  three-foot-depth 
rod  . 


Figure  14  S22  Site  1,  offset  scan  of  three-foot-depth 
rod  . 
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Figure  23  Su  Site  2,  center  scan  of  two-foot-depth 


Figure  24  Su  Site  2,  offset  scan  of  two-foot-depth 
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Figure  27  S22  Site  2,  center  scan  of  two-foot-depth 


Figure  28  S22  Site  2,  offset  scan  of  two-foot-depth 


Figure  29  S21  Site  1,  center  scan  of  two-foot-depth 


Figure  30  S21  Site  1,  offset  scan  of  two-foot-depth 


Figure  35  Su  Site  2,  center  scan  of  one-foot-depth 


Figure  36  Su  Site  2,  offset  scan  of  one-foot-depth 


Figure  41  S21  Site  1,  center  scan  of  one-foot-depth 


Figure  42  S21  Site  1,  offset  scan  of  one-foot-depth 


Figure  47  Su  Site  2,  time-domain  empty  site 


Figure  48  S22  Site  2,  time-domain  empty  site 


Figure  49  S21  Site  1,  time-domain  empty  site 


Figure  50  S21  Site  2,  time-domain  empty  site 
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2.3  Statistical  Characterization  of  Clutter 


As  noted  above,  the  root  response  is  not  distinct  in  the  measurements.  To  better  characterize 
early-time  clutter  in  these  time-domain  GPR  returns,  they  were  examined  statistically.  In  this 
effort  the  41  time-domain  returns  for  the  Sn,  S21,  and  S22  channels  over  the  interval  from  6  ns 
(ground  surface  location)  to  13.5  ns  were  treated  as  realizations  of  stochastic  processes.  The 
standard  deviations  of  these  returns  (also  time  functions)  were  then  computed.  The  implicit 
assumption  here  is  that  scattering  from  randomly  placed  roots  would  tend  to  produce  higher 
standard  deviations  at  depths  where  they  were  present.  The  calculation  was  repeated  for  the 
empty  site,  the  two-foot-depth  rod,  and  the  three-foot-depth  rod.  The  one-foot-depth  data  was 
not  used  because  the  target  response  begins  at  approximately  1 1  ns  and  would  interfere  with  the 
distribution.  The  standard  deviations  of  the  Si  1,  S22,  and  S21  channels  are  shown  in  Figure  51, 
Figure  52,  and  Figure  53,  respectively.  As  noted  previously,  the  surface  scattering  appears  at 
roughly  time=6.5  ns.  The  frequency  range  used  for  the  plots  was  0.05-1.5  GHz.  It  is  interesting 
to  observe  that  the  standard  deviation  values  corresponding  to  the  top  12  inches  of  soil  (i.e.,  prior 
to  t=12  ns)  are  higher  in  magnitude  for  Site  1  (the  high  root  density  site)  than  for  Site  2  (the  low 
root  density  location)  by  6  to  12  dB  for  all  polarizations.  As  noted  above  (and  shown  in  Section 
4)  the  roots  are  found  overwhelmingly  in  the  top  12  inches  of  soil  and,  hence,  the  increased 
clutter  is  evidence  for  root  scattering. 

An  analysis  of  scattering  in  various  frequency  bands  revealed  that  the  largest  difference  in  the 
Site  and  Site  2  scattering  standard  deviations  occurs  between  250  and  750  MHz,  as  shown  in 
Figure  54  through  Figure  56.  At  lower  frequencies  the  dimensions  of  the  roots  are  probably  too 
small  to  cause  significant  scattering.  At  higher  frequencies,  increased  soil  absorption  reduces  the 
penetration  depth,  making  the  radar  insensitive  to  the  roots.  The  standard  deviation  curves  for 
lower  and  higher  frequency  bands  appear  in  Figure  57  through  Figure  62.  It  is  worth  noting  that 
the  largest  percentage  changes  (and  the  weakest  signal  levels)  occur  for  the  cross-polarized 
results. 
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Figure  51  Sn  Standard  deviation  for  the  full  band  (50  to  1500  MHz) 
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Figure  52  S22  standard  deviation  for  the  full  band  (50  to  1500  MHz) 
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Figure  53  S2i  standard  deviation  for  the  frill  band  (50  to  1500  MHz) 
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Figure  54  Su  standard  deviation  for  the  band  250  to  750  MHz 
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Figure  55  S22  standard  deviation  for  the  band  250  to  750  MHz 
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Figure  56  S21  standard  deviation  for  the  band  250  to  750  MHz 
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Figure  57  Su  standard  deviation  for  the  band  50  to  250  MHz 
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Figure  58  S22  standard  deviation  for  the  band  50  to  250  MHz 
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Figure  59  S2i  standard  deviation  for  the  band  50  to  250  MHz 
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Figure  60  Su  standard  deviation  for  the  band  750  to  1500  MHz 
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Figure  61  S22  standard  deviation  for  the  band  750  to  1500  MHz 
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Figure  62  S2i  standard  deviation  for  the  band  750  to  1500  MHz 


3.0  Surface  Scattering  Measurements  from  Elevated  Horn 


Near-surface  roots  could  produce  surface  clutter  that  would  be  visible  to  a  radar  using  an 
elevated  antenna.  For  a  radar  that  views  the  ground  at  an  oblique  angle,  surface  clutter  cannot  be 
separated  from  clutter  at  depth,  which  makes  it  particularly  troubling  and  a  possible  source  of 
confusion.  We  performed  experiments  to  characterize  surface  clutter  at  the  Eglin  site  prior  to  the 
ground-contacting  measurements  described  in  Section  2.  Because  of  the  logistical  difficulties  of 
fielding  an  elevated  radar,  the  scope  of  this  effort  was  limited.  In  this  section  we  review  the 
procedure  used  in  those  measurements  and  the  results  obtained. 

3. 1  Measurement  Setup 

A  commercial  RF  vector  network  analyzer  (HP8753C)  and  a  standard  horn  antenna  (AEL  model 
H-1734)  were  used  for  collecting  reflection  (Sn)  data  over  the  frequency  band  from  0.05  to  2.0 
GHz.  Both  parallel  and  perpendicular  polarization  data  were  collected,  but  no  cross-polarized 
returns  were  acquired.  The  antenna  was  attached  to  the  end  of  a  wood  4x4  beam  and  placed  on 
the  roof  of  a  rental  truck,  as  shown  in  Figure  63.  The  relevant  dimensions  are  shown  in  Figure 
64.  The  radar  control  and  data  collection  was  done  from  a  laptop  computer  located  inside  the 
truck.  Stepped- frequency  data  were  collected  from  0.05  to  2.0  GHz  with  401  steps  and  a  step 
size  of  4.875  MHz.  The  IF  bandwidth  was  set  to  300  Hz,  and  the  system  power  level  was  set  to 
0  dBm. 


Figure  63  Surface  scattering  measurement  setup  with  calibration  sphere. 
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Figure  64  Dimensions  relevant  to  the  surface  scattering  measurement. 


These  data  were  calibrated  using  a  4- inch  diameter  conducting  sphere  placed  on  a  9-inch  high 
Styrofoam  stand  and  positioned  on  the  horn’s  boresight.  In  addition  to  this  calibration 
measurement,  a  “background”  measurement  was  also  made  with  the  horn  pointed  to  the  sky. 

The  latter  measurement  provides  an  estimate  of  the  system  response,  which  was  used  in 
subsequent  processing. 

At  both  Site  1  and  Site  2  four  frequency- swept  measurements  were  performed:  the  parallel- 
polarization  response  of  the  empty  site,  the  perpendicular-polarization  response  of  the  empty  site, 
the  parallel-polarization  response  of  a  calibration  sphere,  and  the  perpendicular-polarization 
response  of  a  calibration  sphere. 

3.2  Results 

The  surface  scattering  data  were  processed  to  remove  the  system  response  by  subtracting  the 
background  (sky)  measurements.  The  data  were  then  transformed  to  the  time-domain.  Time  t=0 
corresponds  to  the  horn  input  connector.  From  Figure  64  we  see  that  the  first  reflection  from  the 
ground  appears  at  t=27.9  ns.  Returns  prior  to  this  time  may  be  contaminated  by  changes  in  the 
background,  which  included  cable  motion  occurring  when  the  operator  got  in  and  out  of  the  truck 
to  remove  the  calibration  sphere.  A  linearly  decreasing  attenuation  (from  -40  dB  to  0  dB)  was 
applied  to  the  first  20  ns  of  the  results  shown  in  the  following  plots  to  emphasize  the  surface 
scattering  phenomena. 

The  background-subtracted  sphere  returns  (i.e.,  sphere  -  no  sphere  data)  are  shown  in  Figure  65 
and  Figure  66  for  parallel  and  perpendicular  polarization,  respectively.  The  sphere  response  is 
the  peak  observed  between  32  and  33  ns.  The  empty  site,  i.e.  no  sphere,  measurements  for  Site  1 
and  Site  2  are  shown  in  Figure  67  and  Figure  68. 

The  results  do  not  show  a  clear  difference  between  the  sites,  and  we  believe  that  they  are 
generally  inconclusive.  There  is  a  large  difference  in  the  early-time  perpendicular  polarization 
return  of  Figure  66,  which  we  attributed  to  the  aforementioned  cable  motion.  The  remaining 
results  do  not  show  a  clear  difference  between  the  two  sites.  Between  45  and  55  ns  the 
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magnitude  of  the  reflected  signal  from  Site  1  is  approximately  5  dB  above  the  signal  from  Site  2, 
but  that  factor  is  not  consistent  over  a  larger  time  interval.  A  larger  sample  size  will  be  required 
before  a  definitive  assessment  of  surface  scattering  can  be  done. 
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Figure  65.  Background-suppressed  data  (i.e.,  sphere  -  no  sphere)  for  the  parallel  polarization  plotted  in  the  time- 

domain  for  the  frequency  band  50  -  2000  MHz. 
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Figure  66.  Background-suppressed  data  (i.e.,  sphere  -  no  sphere)  for  the  perpendicular  polarization  plotted  in  the 

time-domain  for  the  frequency  band  50  -  2000  MHz 


X 


MsgnitHte-K®!  MflgnHudeHdB) 


PufilM  Pi.-lSriidlitn  (EirluhfJ 


-50 - - - t - : - : - r -  -. - T - =^= - 

-  Site  f 


0  ID  20  30  4Q  50  00  70  &Q  SO  TOO 


True  fns| 

Figure  67  Parallel  polarization  return  for  the  empty  site  (50  -  2000  MHz) 
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Figure  68  Perpendicular  polarization  return  for  the  empty  site  (50  -  2000  MHz) 


4.0  Supplemental  Measurements 


In  addition  to  the  scattering  measurements  documented  in  Sections  2  and  3,  measurements  of 
various  site  characteristics  were  made.  In  this  section  we  document  those  measurements. 

4. 1  Root  Structure  and  Properties 

After  completing  all  radar  measurements,  Sites  1  and  2  were  excavated  to  a  depth  of  24  inches. 
Care  was  taken  to  preserve  any  roots  found.  Figure  69  shows  the  root  structures  found  at  Sites  1 
and  2.  In  both  cases  the  roots  have  been  painted  red  to  increase  their  photographic  contrast.  It  is 
evident  that  there  are  essentially  no  roots  present  at  Site  2.  Since  the  sample  locations  had  to  be 
chosen  prior  to  their  excavation,  this  finding  was  a  welcome  discovery.  It  also  suggests  that 
surface  shrub  density  is  a  good  predictor  of  subsurface  density. 


(a)  Site  1  (b)  Site  2 

Figure  69.  The  test  sites  after  excavation  to  a  depth  of  24  inches.  The  roots  were  left  intact  and 
painted  red  for  better  visibility.  String  lines  mark  1  meter  square  sections 

Several  features  of  the  roots  at  Site  1  require  additional  discussion.  A  number  of  overlapping 
photos  of  that  location  were  taken,  and  the  results  were  ortho-rectified  to  produce  the  plan- view 
image  shown  in  Figure  70,  in  which  the  image  distances  approximate  the  true  plan  distances. 

The  large  root  mass  shown  near  the  center  of  the  figure  was  not  evident  from  the  surface,  but  it 
was  probably  the  nucleus  of  a  vigorous  shrub.  Roots  radiate  from  it  in  essentially  all  directions. 
That  object,  which  is  approximately  10  cm  in  diameter,  comprises  the  largest  object  seen  in  our 
study,  and  it  features  prominently  in  the  analysis  in  Part  2  of  this  report.  Unfortunately,  as 
indicated  in  the  figure,  the  GPR  scan  did  not  traverse  the  large  root  mass.  (The  “centered”  scan 
path  is  shown.  The  “offset”  path  would  be  below  the  bottom  edge  of  the  image.) 

The  entire  root  structure  shown  in  Figure  70  was  exhumed  and  returned  to  our  Laboratory  for 
additional  study.  Root  diameters  were  measured,  and  a  histogram  of  the  results  (normalized  per 
unit  length)  appears  in  Figure  71.  Very  small  roots  were  common,  but  did  not  survive  excavation 
and  are  unlikely  to  affect  the  GPR  measurements.  The  distribution  suggests  a  modal  diameter  of 
about  1.5  cm.  Several  root  samples  were  also  excavated  and  immediately  sealed  for  further 
analysis. 
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Figure  70.  Top  view  of  the  root  distribution  at  Site  1.  The  approximate  radar  path  is  shown  by  the 

dashed  line. 
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Figure  71  Histogram  of  root  diameters. 


Soil  samples  were  obtained  every  six  inches  to  a  depth  of  three  feet.  Those  samples  were 
double-bagged  in  sealed  plastic  bags  and  returned  to  our  facility  for  moisture  measurements  and 
dielectric  property  measurements.  Results  of  those  measurements  appear  in  Part  2  of  this  report. 
The  local  soil  was  very  sandy,  but  the  top  surface  had  a  somewhat  different  composition.  Figure 
72  shows  that  roughly  the  top  six  inches  of  soil  has  a  markedly  darker  color  than  the  deeper  soil, 
possibly  as  a  result  of  decayed  organic  matter. 


Figure  72.  Near-surface  soil  shows  darker  color  than  deeper  soil. 


4.2  Shrub  density. 

Finally,  the  spatial  density  of  shrubs  in  the  region  of  our  experiments  was  measured.  A  total  of 
six  30  foot  square  regions  (arranged  as  a  2x3  grid)  were  marked,  and  the  shrubs  within  those 
areas  were  counted.  The  counts  appear  in  Table  2,  which  shows  a  mean  density  of  one  shrub  per 
30.5  square  feet,  or  a  mean  distance  of  5.5  feet  between  shrubs. 


xiv 


Table  2.  Shrub  counts  for  adjacent  30  foot  square  regions. 


Region 

Count 

1 

25 

2 

21 

3 

40 

4 

26 

5 

39 

6 

26 
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5.0  Summary  and  Conclusions 


A  combination  of  experiment  and  theory  was  used  to  study  scattering  from  buried  roots.  In  this 
document,  we  present  the  experimental  aspects  of  the  work.  A  number  of  measurements  were 
performed  by  the  Ohio  State  University  ElectroScience  Laboratory  (OSU  ESL)  at  Site  C62, 

Eglin  AFB,  FL.  Wideband,  multi-polarization  GPR  measurements  were  made  with  a  radar  that 
uses  a  ground-contacting  antenna.  The  frequency  band  examined  in  this  work  was  50-2000 
MHz,  although  unexpected  clutter  in  one  polarization  forced  us  to  use  only  results  below  1500 
MHz. 

Some  interesting  results  were  gleaned  from  time-domain  (“waterfall”)  plots  of  the  radar  output 
as  it  is  scanned  over  sites  with  and  without  roots.  A  comparison  of  the  co-polarized  returns  for 
sites  with  and  without  roots  showed  very  similar  phenomena,  including  scattering  from  largely 
random  objects  in  the  top  foot  of  soil.  These  scattering  events  were  weaker  than  the  return  from  a 
steel  rod  inserted  into  the  site  from  a  distant  access  trench.  A  statistical  analysis  of  the  returns 
showed  that  the  site  with  roots  had  a  generally  higher  clutter  level  than  the  site  without  roots. 

The  difference  between  the  sites  is  perhaps  most  evident  for  the  cross-polarized  returns.  Sites 
with  roots  produce  a  significantly  more  cross-polarized  return  than  sites  without  roots.  The 
statistical  analysis  confirmed  this  finding.  For  both  co-polarized  and  cross-polarized  returns,  root 
clutter  was  found  to  be  strongest  in  the  frequency  band  from  250  to  750  MHz. 

Measurements  of  surface  clutter  were  also  made  using  an  elevated  antenna.  Because  of  the 
difficulty  of  fielding  an  elevated  radar,  the  scope  of  those  measurements  was  limited.  An 
analysis  of  the  results  yielded  no  obvious  difference,  but  it  is  difficult  to  draw  compelling 
conclusions  from  them  because  of  the  small  data  set  (only  two  sites  were  examined.) 

Additional  investigation  of  this  scattering  from  biomass  appears  in  Part  2  of  this  report.  In  that 
document  we  describe  a  numerical  model  developed  by  OSU  and  capable  of  predicting  many  of 
the  results  shown  herein. 
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Appendix:  File  Formats 

GPR  data  files 


The  data  for  system-calibrated  frequency-domain  radar  data  was  stored  in  binary  files  that  use 
the  naming  convention  “aydddaa  .  cdt”,  where  “a”  is  a  letter  in  the  range  a-z  used  for  file 
ordering.  The  letter  “y”  indicates  the  last  digit  of  the  year.  For  example,  “0”  represents  the  year 
of  2000.  The  three-digit  number,  “ddd”  ,  indicates  the  Julian  date  when  the  data  was  stored. 
Each  file  contains  two-dimensional  (space  and  frequency)  data  taken  at  multiple  target  locations. 
Frequency  responses  are  stored  column-wise  with  the  first  frequency  (10  MHz)  stored  in  the  first 
row,  the  second  frequency  (12  MHz)  data  stored  in  the  second  row,  etc.  Table  3  illustrates  the 
new  format  of  multiple  GPR  data.  Data  collected  at  each  position  is  contained  in  one  vertical 
block  that  comprises  six  columns  corresponding  to  the  real  and  imaginary  parts  of  Sn,  S21  and 
S22,  respectively,  as  illustrated  in  Table  4.  For  each  radar  file,  "*  .  cdt”,  there  is  an  associated 
comment  file  called  “*  .  txt”  to  store  the  system  information  and  comments.  There  is  also  an 
associated  processed  file  called  “*.mat”  file  that  stores  the  processed  results  and  all  processing 
parameters.  All  of  these  files  are  available  upon  request. 


Table  3  Format  for  multiple-position  GPR  data. 


Data(fi,Xi) 

Data(fi,X2) 

Data(fi,Xn) 

Data(f2,Xi) 

Data(f2,X2) 

Data(f2,Xn) 

Data(fm,Xi) 

Data(fm,X2) 

Data(fm,Xn) 

An  example  of  the  format  for  Data(fm,Xn)  is  shown  below. 

Table  4  Data  format  for  a  single  frequency  at  a  single  position. 


Re(Sn) 

Im(Sn) 

Re(S21) 

Im(S2i) 

Re(S22) 

Im(S22) 

"Re()"  and  "Im()"  indicate  the  real  and  imaginary  parts,  respectively. 

Comment  (Text)  Files 

The  comment  text  files  contain  information  about  measurement  conditions  (i.e.  position, 
direction,  etc.)  and  any  comments  the  user  entered  during  the  measurements.  For  example,  the 
comment  file  a0348ag.txt  is  printed  below: 

Start  Time:  511/Stop  Time:  519/Target  #:  A2/Start  f:  10/Stop  f:  410/Num  of 
Points:  201/Antenna  Orientation:  92/Num  of  Scans:  41/Antenna  Position: 
xstart:  -60/delta  x:  3/yoffset:  0/Relative  Permitivity:  20/User  Comments:  /$ 

The  file  contains  the  time,  in  minutes  after  midnight,  at  which  the  scan  was  started  and  finished. 
It  also  contains  the  information  necessary  to  recreate  the  frequency  vector  and  the  position 
vector,  i.e.  number  of  positions  of  data  taken  and  the  spacing  between  them. 
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1.0  Introduction 


This  document  comprises  Part  2  of  the  final  report  for  the  research  effort  entitled 
“GPR  performance  in  the  presence  of  buried  biomass”  (Contract  DAAD 19-01- 1-0760, 
Project  No.:  E-43042-EL-000-0 1234-1).  There  is  extensive  anecdotal  evidence  that 
biomass  (i.e.,  plant  and  tree  roots)  adversely  affects  ground-penetrating  radar  (GPR) 
performance.  In  this  effort  the  Ohio  State  University  (OSU)  ElectroScience  Laboratory 
(ESL)  described  the  physics  that  govern  the  propagation  and  scattering  of  GPR 
waveforms  in  the  presence  of  buried  biomass.  This  work  is  being  performed  in 
collaboration  with  the  US  Army  Research  Laboratory. 

This  report  is  presented  in  two  parts.  In  Part  1  we  document  the  results  of 
experimental  GPR  measurements  performed  at  Site  C62  at  Eglin,  AFB,  FL.  In  Part  2  we 
describe  numerical  modeling  work  that  permits  us  to  simulate  and  interpret  many  of  the 
results  presented  herein. 

Ground  penetrating  radar  (GPR)  has  been  used  extensively  to  detect  and 
discriminate  unexploded  ordnance  (UXO),  because  of  its  ability  to  penetrate  soil.  Recent 
data  suggest  that  the  performance  of  GPR  may  be  impaired  by  buried  biomass. 
Specifically,  problems  were  observed  during  measurements  by  the  U.S.  Army  Research 
Laboratory  (ARL)  at  Eglin  AFB  using  the  BoomSAR~an  ultra-wideband  (UWB) 
synthetic  aperture  radar  (SAR)  mounted  on  a  mobile  lift  platform.  It  is  known  that  there 
is  a  high  density  of  roots  at  that  site,  and  those  roots  were  thought  to  affect  the 
penetration  of  electromagnetic  energy  into  the  soil. 

A  team  from  OSU/ESL  visited  Eglin  AFB  during  December  10-12,  2001  to 
perform  additional  tests.  A  detailed  discussion  of  the  OSU  measurements  appears  in  a 
separate  report  [1].  GPR  briefly,  measurements  were  performed  in  areas  roughly  2  m  by  2 
m  in  size.  One  of  those  areas  had  a  high  density  of  roots  and  is  referred  to  as  “Site  1.” 
This  area  is  shown  in  Figure  1.1  after  the  soil  was  removed.  Most  of  the  roots  were 
found  in  the  top  12  inches  of  soil.  The  measurements  comprised  scans  over  the  roots  with 
a  ground-contacting  GPR  antenna. 
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The  root  structure  shown  in  Figure  1.1  was  removed  and  taken  to  OSU  for  further 
study.  All  of  the  roots  simulations  presented  here  are  based  on  the  root  structure  shown  in 
that  figure. 


Figure  1.1.  “Site  1”  area 

In  this  report  we  present  a  combination  of  experiments,  theory  and  modeling  that 
was  used  to  gain  additional  insight  into  root  scattering.  The  report  is  organized  in  five 
major  chapters.  In  Chapter  2  we  describe  the  discrete  dipole  approximation  (DDA),  an 
algorithm  that  can  be  used  to  compute  fields  scattered  from  tree  roots.  In  Chapter  3  the 
exact  solutions  and  small  argument  approximations  for  scattering  from  infinitely  long 
cylinders  are  reviewed.  The  use  of  the  DDA  as  a  model  for  the  GPR  system  is  presented 
in  Chapter  4.  Case  studies  of  the  effects  of  roots  on  GPR  are  presented  in  Chapter  5.  The 
report  concludes  with  a  discussion  of  the  results  in  Chapter  6. 
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2.0  Computing  Root  Scattering 


Scattering  from  buried  objects  has  been  studied  extensively,  because  of  its  many 
applications  including  biomedical  imaging,  remote  sensing,  and  geological  exploration 
[2].  The  targets  of  interest  are  both  dielectric  and  conducting  objects.  Many  papers 
relevant  to  this  topic  have  been  published  in  the  last  few  decades,  and  only  a  few  key 
works  will  be  cited  here.  Homogeneous  dielectric  cylinders  in  a  lossy  half-space  were 
considered  by  Mahmoud  et  al.  [3] .  Zhuck  et  al.  [4]  investigated  inhomogeneous  dielectric 
cylinders  embedded  in  a  stratified  medium.  Three-dimensional  (3-D)  objects  buried  in  a 
dielectric  half-space  were  studied  by  Lombardi  et  al.  [2].  Butler  et  al.  investigated 
conducting  cylinders  of  arbitrary  cross  sections  partially  buried  between  two  media 
[5], [6].  Perfectly  conducting  objects  of  arbitrary  shape  in  multilayered  media  were 
explored  by  Michalski  et  al.  [7].  3-D  dielectric  and  conducting  objects  buried  under  a 
one-dimensional  (1-D)  multilayered  medium  were  studied  by  Cui  et  al.  [8].  Since  buried- 
targct  problems  involve  Sommerfeld-like  integrals,  which  are  time-consuming  to 
evaluate,  fast  algorithms  have  been  developed,  e.g.,  Cui  et  al.  [9], [10]. 

Although  many  studies  of  scattering  from  dielectric  objects  under  a  half-space 
have  appeared,  few  of  these  specifically  address  tree  roots.  Sullivan  [11]  modeled  tree 
roots  as  an  ensemble  of  connected  cylinders  without  including  interaction  between  each 
root.  El-Shenawee  et  al.  [11]  modeled  a  dielectric  clutter  object  (which  could  be  a  tree 
root)  close  to  an  antipersonnel  mine-like  object  under  a  rough  surface. 

In  the  present  chapter  we  present  a  novel  approach  to  the  root  scattering  problem. 
A  formulation  of  the  underlying  integral  equation  and  its  discretization  based  on  the 
“Discrete  Dipole  Approximation  (DDA)”  [13]  are  described.  Computational  details  are 
also  discussed.  A  validation  of  the  DDA  code  is  presented,  and  its  limitations  are 
explained. 

2.1  Formulation  of  Volume  Integral  Equation 

As  noted  in  Chapter  1,  all  of  the  analysis  presented  herein  is  based  on  the  root 
structure  shown  in  Figure  1.1.  The  root  structure  has  been  modeled  by  an  ensemble  of 
cylindrical  segments.  A  comparison  of  the  original  and  discretized  versions  is  shown  in 
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Figure  2.1.  The  cylinder  diameters,  which  range  from  0.3  -  4  cm,  were  obtained  from 
measurements  of  the  structure.  A  large  root  mass  appears  near  the  center  of  the  figure. 
For  reasons  that  will  be  discussed  later,  that  mass  was  represented  by  a  group  of  densely 
packed  cylinders. 
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Figure  2.1.  Comparison  between  the  original  and  discretized  root  structures,  (a) 
Experimental  image,  (b)  Model  representation. 

Each  root  segment  is  described  by  several  variables,  and  from  Figure  2.1(b)  it  is 
evident  that  we  have  to  deal  with  a  large  number  of  unknowns.  Also,  we  try  to  simulate 
GPR  data,  which  require  computation  of  fields  at  multiple  frequencies  and  positions.  This 
led  us  to  investigate  efficient  algorithms.  Note  that  herein,  we  will  neglect  the  air-ground 
interface  (explain  in  Chapter  4)  so  that  the  CPU  and  memory  requirements  can  be 
reduced. 

In  the  following  a  formulation  of  the  electric-field  volume  integral  equation  is 
discussed  using  a  discretization  based  on  the  DDA  [13], [14].  The  DDA  is  used  because 
of  its  efficiency  in  computing  fields  scattered  from  complicated  dielectric  bodies. 
Furthermore,  it  includes  interactions  all  parts  of  the  bodies.  In  using  the  DDA  we  will 
approximate  a  continuum  dielectric  target  by  a  finite  array  of  dipoles.  Hence,  we  expect 
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accurate  results  only  if  the  array  can  describe  the  shape  of  the  target  satisfactorily.  Tests 
are  presented  that  quantify  the  algorithm’s  performance. 


V m>  \£m>fJ‘o) 


Figure  2.2.  A  plane  wave  incident  on  a  scatterer  with  permittivity  £,  permeability  ju0  and 
volume  V.  The  scatterer  resides  in  a  region  denoted  Vm,  which  contains  a  medium  with 
parameter  ( £m  ,  jU0). 


Consider  a  plane  wave  incident  on  a  scatterer  with  permittivity  £  (Figure  2.2).  For 
the  above  geometry,  Maxwell’s  equations  can  be  written  as 


VxE  =  +ia/ioH 

(2.1) 

VxH  =-ioj£  p(r)E 

(2.2) 

where 


ep(r)  = 


£ 


if  r  e  V 
if  re  Vm  ' 


(2.3) 


The  time  dependence  e',ax  is  implied  and  suppressed  through  out  the  derivation. 

Taking  the  curl  of  equation  (2.1),  using  equation  (2.2)  and  subtracting  k2E  yields  the 
vector  Helmholtz  equation 

VxVxE-k2E=k2(£c -\)E  (2.4) 


where  £c  =  £  I  £m  and  k 2  =  of  /Llo£m  .  Throughout  this  work,  the  term  “relative  permit¬ 
tivity”  will  refer  to  £c ,  which  is  computed  for  a  background  medium  £m .  The  well- 
known  solution  of  equation  (2.4)  can  be  expressed  as 

E{f)  =Emc  +J  d3xG(r,r)  k2(£c  -1  )E(r')  (2.5) 
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where 


G(r,r')  =  (/  +  ^)  g(r,r')  (2.6) 

k 

and 

—/A:  |  r—r'  | 

g(r,r')  =  f_  _  (2.7) 

r  —  r 

The  dyadic  Green’s  function  G  has  a  singularity  given  by  [15], [16] 

—  —  J  fit  y  —  r/'\ 

G{f,r)  =  PVG{f,r) - ^ - -  (2.8) 

k 

where  PV  is  an  abbreviation  for  the  principal  value,  and  S(r-r')  is  the  three- 
dimensional  Dirac  delta  function.  The  integral  of  the  first  term  in  (2.8)  is  performed  over 
the  volume  of  the  scatterer  excluding  an  infinitestimal  volume  Vs  about  the  observation 
point  r  .  Although  V^is  infinitestimal,  its  shape  must  be  specified,  and  that  shape  affects  the 

dyad  L  .  Substituting  (2.8)  into  (2.5),  the  volume  integral  equation  becomes 

E(r)  =  Einc(r)  +  f  d3x  <f(r,r')  •  (£c -1)  k2  E(r)  -  (£c -1)1  ■  E(r) .  (2.9) 

v-vt 

Equation  (2.9)  will  now  be  discretized  and  evaluated  numerically  using  the  DDA 
algorithm. 

The  induced  polarization  density  P  can  be  expressed  as  a  function  of  E  as 

P(r)  - 

-A-L  =  (sc-l  )E(r).  (2.10) 

Gn 

The  volume  V  is  discretized  into  volumes  AVj  with  centers  at  ry ,  j  =1,  2,...,  N  (N  is  the 
number  of  dipoles).  The  field  Ej  inside  AVj  is  further  approximated  as  a  constant  (“pulse 
expansion”)  to  obtain 

P(r)  Pj  - 

^1  =  -l  =  {£c-\ )E(r),  re  AVj.  (2.11) 

£m  £m 

Then, 
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E(r, )  =  E"*(r, )  +  —  £  G(r.Jf)-  P.\VJ  + 


P 

C/m  j=hj*i 


~(ee-l)  LEifi) 


k2  f 
vk 

(2.12) 


where  we  have  approximated  the  PV  integral  away  from  r  by  its  mean  value,  and  we 

have  taken  PV |()<i3v  ~  |()<i3.x:  for  a  finite  Vs-  Let  R  be  the  unit  vector  from  f;  to  r , 

Vi  v-vs 


let/?  =|  r;  -ry|  and  let  I  be  the  unit  dyad.  For  ri  rj9  G  ( ,  r) )  can  be  found  from  (2.6) 
by  direct  differentiation 


(ri,rj)  =  {(- 


ikR 


1  +  ikR  +  k2R2)I  +  (3-3  ikR-kzRz)RR 


2  2  \ 


2  r>3 


4  nk-R 


(2.13) 


Define  A(r,r)  =  — —  G(r,r),  D  =  f  d3xG(r.,r') ,  and  let  the  polarization  p.  be 

V-Vs 

given  by  P.AV. .  Using  the  relation  in  (2.1 1),  (2.12)  becomes 


CC~Vpi+  ftA(rl,rJ)-pJ  =  E»!(rl)  (2.14) 

j=hj*i 


where  a 


1 

£m(ec  -DAV, 


f-  p 

1 

I-k2 

D - - 

l  k  J 

Vk 

1 

1 _ 

If  A  (r^r^)  is  defined  as  ff 


(2.14)  can  be  written  as 

-  Emc  (r) .  (2.15) 

j= i 

From  (2.15)  it  is  evident  that  the  scattering  problem  reduces  to  finding  the 
polarizations  pj  that  satisfy  a  system  of  3N  complex  linear  equations  [A]  [p]  =  ( E ""  J . 
Once  pj  is  found,  one  can  approximate  the  scattered  fields  by 


El(r)»ftA(r,rj)-pJ.  (2.16) 

j= i 

The  polarization  p j  in  (2.15)  can  be  found  only  if  D  and  L  are  known.  Note  that 
both  D  and  L  depend  on  the  shape  of  the  exclusion  volume.  Since  scattering  from 
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cylinders  is  of  interest  here,  it  is  natural  to  take  the  shape  of  the  exclusion  volume  to  be  a 
circular  cylinder. 


Figure  2.3.  (a)  A  thin  cylinder  with  radius  a  and  length  L.  (b)  Segmentation  of  the 
cylinder  into  small  cells  of  length  /. 


Calculation  of  these  dyads  proceeds  as  follows:  First,  consider  a  thin  cylinder 
with  radius  a  and  length  L  »  a  as  shown  in  Figure  2.3(a).  If  the  cylinder  is  also 
electrically  thin  with  (ka)  «  1 ,  the  current  over  each  cross  section  is  approximately 
uniform,  and  only  the  longitudinal  variation  in  pis  of  interest.  The  cylinder  may  be 
segmented  into  small  cells  of  length  l  as  shown  in  Figure  2.3(b).  We  define  an  exclusion 
volume  Vs  to  be  a  cell  with  radius  a  and  length  /.  It  can  be  shown  [14]  that  D  and  L  for 


a  circularly  cylindrical  cell  with  axis  z  are  as  follows. 

L  =  L:xx  +  Ly  yy  +  Lz  zz 

(2.17) 

D  =  Dxxx  +  Dy  yy  +  D.zz 

(2.18) 

L  -L  ~  ^ _ 

'  2a/4  a2 +l2 

(2.19) 

Z  v/4  a2+l2 

(2.20) 

D.  =D.  =^^  +  -{jl2  +4a2  -/}+— In 


ikccl  l 

- 1 - 

6  8 


1  +  -J  l~  +  4cr 
2  a 


D  =  6 


ikal  a2  ,  f  l  +  y]l2  +  4r/ "  ^ 


+  — In 

2 


V 


2  a 


(2.21) 

(2.22) 
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2.2  Computation  of  A  (rt ,  r. ) 

For  ij  ^  r.  the  above  formulation  assumes  that  the  matrix  element  are  given  by 

Au=—  \d3x'G(ri,7)»  —  G(fi,rj)AVj  =  A(ij,r,).  (2.23) 

^ m  A v  j  Gh 

This  approximation  is  not  accurate  if  n  and  r.  are  within  some  small  distance  r,  [  1 7J . 
One  can  improve  the  accuracy  by  applying  numerical  integration  over  the  volume  V) 
centered  at  r. .  An  improved  evaluation  of  the  matrix  terms  are  then  given  by 


\drA(rnr) 

<  j  Vj 

for  r;  -  ry.  \  <rd 

(2.24) 

C  J 

for\r-rj  \>  rd 

For  a  thin  cylinder,  a  suitable  limit  is  rd  =  Xm  /  8  where  Xm  is  the  wavelength  in  the 
medium. 

Figure  2.1(b)  indicates  that  many  cylinders  are  necessary  to  represent  the  root 
structure,  which  corresponds  to  a  very  large  number  of  dipoles.  For  N  dipoles,  storage  of 

the  complex  matrix  [A  ]  requires  3./V x 3iV x 8 bytes  for  four-byte  floating  point  values. 
As  an  example,  2000  dipoles  requires  288  MB  of  memory. 

To  reduce  the  memory  usage,  only  those  terms  A  /y  that  require  numerical 

integration  are  stored.  Terms  that  can  be  expressed  in  closed  form  are  recalculated  every 
time  they  are  used. 


2.3  Solving  Equation  (2.15) 

The  system  of  complex  linear  equations  in  equation  (2.15)  must  be  solved  to  find 
the  polarization  py .  An  iterative  solver  is  used  here  instead  of  a  direct  method,  because 


iteration  is  more  efficient  when  the  number  of  dipoles  is  large.  The  Bi-Conjugate 
Gradient  Stabilized  (Bi-CGSTAB)  method  was  chosen  here,  since  it  has  been  shown  to 
converge  more  rapidly  than  other  methods  (for  DDA)  [18].  The  iterative  solver 


terminates  when  achieving 


[A]  ■  [p]~\Einc] 


<  e  Emc  ,  where  e  is  an  error  criterion, 


chosen  to  be  10"  in  this  work. 
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2.4  Representing  Thick  Cylinders  with  DDA 

The  thin-cylinder  criterion  used  in  DDA  requires  that  the  length  and  diameter  of 
each  cell  must  be  small  compared  to  the  wavelength  of  the  incident  field.  Thus,  we 
require 

^\k(2anrJ<l  (2.25) 

and 

a/7~  kl  <1.  (2.26) 

For  a  cylinder  with  radius  a  >  amax,  a  number  of  thin  cylinders  of  the  same  length 
can  be  used  to  approximate  the  larger  one  such  that  the  volume  is  still  the  same. 
Examples  of  how  to  arrange  the  component  cylinders  are  shown  in  Figure  2.4. 


-0.04  -0.03  -0.02  -0.01  0  0.01  0.02  0.03  0.04  -0.04  -0.03  -0.02  -0.01  0  0.01  0.02  0.03  0.04 


(a)  (b) 

Figure  2.4.  Approximations  for  cylinders  with  different  radii  for  amax  =  0.015  m.  (a)  a  = 
0.030  m.  (b)  a  =  0.0396  m. 

2.5  Validation  with  Other  Codes 

The  DDA  code  has  been  compared  with  the  Boeing  “Body  of  Revolution  (BOR)” 
code  for  cylinders  with  different  radii.  That  code  is  based  on  research  by  Putnam  et  al. 
[19],  Andreasen  [20],  and  Mautz  et  al.  [21].  The  problem  geometry,  shown  in  Figure  2.5, 
comprises  a  plane  wave  impinging  on  a  cylinder  at  different  angles  for  vertical  and 
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horizontal  polarizations.  Figures  2.6,  2.7  and  2.8  show  good  agreement  of  co-polarized 
(VV  and  HH)  backscattered  radar  cross  sections  (RCS)  for  the  DDA  and  BOR  codes, 
even  when  thick  cylinders  are  approximated  by  an  ensemble  of  thinner  cylinders. 


Figure  2.5.  Plane  waves  incident  on  a  cylinder  with  angle  0. 


(a)  (b) 

Figure  2.6.  Comparisons  of  radar  cross  sections  computed  by  DDA  and  BOR  at 
frequency  1  GHz  for  a  0.3  m  long,  0.015  m  radius,  cylinder  with  permittivity  2  (1 
component  cylinder,  i.e.,  a  <  amax).  (a)  VV.  (b)  HH. 
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VV  Radar  Cross  Section  (dB-m2)  VV  Radar  Cross  Section  (dB-m2) 


9  (degree)  9  (degree) 


(a)  (b) 

Figure  2.7.  Same  as  Figure  2.6  but  with  radius  =  0.030  m  (4  component  cylinder,  see 
Figure  2.4(a)).  (a)  VV.  (b)  HH. 


0  (degree)  ®  (degree) 

(a)  (b) 

Figure  2.8.  Same  as  Figure  2.6  but  with  radius  =  0.0396  m  (7  component  cylinder,  see 
Figure  2.4(b)).  (a)  VV.  (b)  HH. 


To  quantify  the  efficiency  of  the  DDA,  we  compared  computational  time  with 
other  codes.  The  geometry  of  the  problem  is  exactly  the  same  as  Figure  2.5  with  0  =  90 
degrees.  The  incident  field  is  vertically  polarized.  The  dimension  of  the  cylinder  is  0.2  m 
in  length  and  0.0667  m  in  radius  with  a  dielectric  constant  of  2.  We  computed  the 
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frequency  response  of  this  cylinder  from  0.05  to  1.50  GHz.  Since  the  BOR  code 
mentioned  before  does  not  scan  in  frequency,  the  OSU’s  Electromagnetic  Surface  Patch 
(ESP)  [22]  code  was  used  instead.  Figure  2.9  shows  good  agreement  of  backscattered 
radar  cross  sections  computed  by  the  DDA  and  ESP  codes.  The  time  it  takes  for  them  to 
compute  the  above  results  is  17  minutes  and  7.9  hours,  respectively,  on  a  Pentium  4  1.60 
GHz  computer  with  256  MB  memory.  Note  that  this  is  probably  not  a  fair  comparison, 
because  DDA  uses  an  iterative  method  to  solve  the  linear  equations  while  ESP  employs 
LU  decomposition  of  the  impedance  matrix.  Nonetheless,  the  large  difference  in  time 
suggests  that  DDA  is  an  obvious  choice  for  handling  large  problems. 


Figure  2.9.  Comparison  of  VV  radar  cross  section  computed  by  DDA  and  ESP.  The 
geometry  of  the  problem  in  shown  in  Figure  2.5  with  0  =  90  degrees.  The  cylinder  is  0.2 
m  in  length  and  0.0667  m  in  radius. 

2.6  Limitations  of  DDA 

DDA  has  a  few  limitations,  which  are  now  described.  First,  in  its  current  form 
DDA  can  only  be  used  to  compute  fields  scattered  from  a  dielectric  target  with  a 
dielectric  constant  that  is  not  too  large  [13].  It  has  been  found  that  the  code  is  accurate 
when  the  relative  permittivity  magnitude  is  less  than  5  to  6. 
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Second,  there  should  be  no  intersections  between  the  dielectric  rods.  Such 
intersections  can  result  if  dipoles  from  one  cylinder  are  too  close  to  dipoles  from  another 
cylinder.  To  prevent  this,  a  small  gap  between  adjacent  cylinders  is  necessary,  and  Figure 
2.1(b)  shows  how  small  gaps  are  included  in  the  discretized  root  model.  It  is  reasonable 
to  question  whether  this  gap  produces  significant  errors  in  the  calculation,  we  now  show 
that  these  errors  can  be  neglected. 

The  responses  from  a  2  cm-diameter  cylinder  with  and  without  a  1  cm  gap  are 
calculated.  The  geometries  of  these  two  cases  are  shown  in  Figures  2.10(a)  and  2.10(b), 
respectively.  Note  that  the  background  is  free  space.  Figures  2.11  and  2.12  compare  the 
backscattered  responses  of  the  geometries  in  both  the  frequency  and  time  domains  for  the 
broadside  incident  fields. 


z  z 

A  A 


Figure  2.10.  Cylinder  geometries  used  to  test  the  effects  of  gaps,  (a)  Cylinder  with  a  1  cm 
gap.  (b)  Cylinder  without  a  gap. 
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HH  Radar  Cross  Section  (dB-m2)  VV  Radar  Cross  Section  (dB-m2) 


(a)  (b) 

Figure  2.11.  Comparison  of  fields  backscattered  from  a  cylinder  with  and  without  a  gap 
for  relative  permittivities  of  2  and  4.  The  polarization  of  the  incident  field  is  vertical,  (a) 
Frequency  domain  (RCS).  (b)  Time  domain. 


(a)  (b) 

Figure  2.12.  Same  as  Figure  2.11  but  for  horizontal  polarization,  (a)  Frequency  domain 
(RCS).  (b)  Time  domain. 
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The  highest  frequency  considered  here  produces  a  wavelength  that  is  comparable 
to  that  produced  at  1  GHz  in  soil  with  a  relative  permittivity  of  5.  From  Figures  2.11  and 
2.12,  it  is  clear  that  the  effect  of  a  gap  can  be  neglected. 

2. 7  Summary 

A  root  structure  found  at  Eglin  AFB,  FL,  was  represented  by  an  ensemble  of 
dielectric  cylinders  for  computational  purposes.  The  formulation  of  the  electric-field 
volume  integral  equation  and  its  discretization  based  on  the  discrete  dipole  approximation 
(DDA)  were  described.  The  DDA  code  was  validated  against  a  body  of  revolution  (BOR) 
code  and  OSU’s  Electromagnetic  Surface  Patch  (ESP)  code.  It  was  shown  that  DDA  can 
predict  fields  scattered  from  dielectric  cylinders  accurately.  The  efficiency  of  the  DDA  is 
illustrated  by  comparing  the  computational  time  with  ESP.  Limitations  of  DDA  were 
discussed  as  well. 
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3.0  Scattering  from  Long  Roots  Approximated  by  Infinitely 
Long  Cylinders 


In  the  previous  chapter  it  was  shown  that  the  DDA  code  can  accurately  compute 
scattering  from  finite,  circular,  dielectric  cylinders.  In  this  chapter  we  will  study  the  exact 
solutions  for  scattering  from  a  dielectric  cylinder  of  infinite  length  and  their 
approximations  for  a  thin  cylinder.  The  above  approximations  are  interesting  not  only 
because  tree  roots  usually  resemble  long,  thin  cylinders,  but  also  because  they  show  how 
the  scattered  field  depends  on  the  radius  and  relative  permittivity  of  the  cylinder.  Here, 
the  cylinders  reside  in  a  lossless  homogeneous  medium  and  normally  incident  plane 
waves  with  TMZ  and  TEZ  polarizations  are  assumed.  Comparisons  of  scattered  fields  for 
these  two  cases  will  be  presented. 

3.1  Scattering  from  a  Dielectric  Circular  Cylinder:  TM-  Polarization 

Assume  that  a  TMZ  plane  wave  is  normally  incident  on  a  dielectric  circular 
cylinder  with  radius  a,  permittivity  £,i  and  permeability  fio  as  illustrated  in  Figure  3.1.  The 
time  dependence  e~uot  is  implied  and  suppressed  throughout  the  derivation. 

y 


Einc  © - ► 


Figure  3.1.  TMZ  plane  wave  impinges  on  a  dielectric  circular  cylinder. 

It  can  be  shown  [23]  that  for  this  geometry,  the  expression  for  the  scattered  fields  is 

E^zE^cX^Me-1’*  (3.1) 

n=— oo 
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where 


/  _  / 


Cn  = 

(  „  Jn  (Pmd)J n(Pda)~  JfTcJ n(Pma)J n  ( Pda ) 

■J£jXMHH'(/}ma)-J„(MH?(/3ma) 

(3.2) 

J,X-) 

=  the  Bessel  function  of  the  first  kind  of  order  n  (n  is  an  integer). 

ja-: 

)  =  derivative  of  Jn{.)  with  respect  to  its  argument. 

**?( 

.)  =  the  Hankel  function  of  the  first  kind  of  order  n. 

ha 

(.)=  derivative  of  H^\.)  with  respect  to  its  argument. 

E0 

=  magnitude  of  incident  fields  (set  to  1). 

Pm 

=  WA)  £,n 

(3.3) 

Pd 

=  £d 

(3.4) 

£ 

—  Sd  /  £m 

(3.5) 

As  noted  previously,  the  term  “relative  permittivity”  will  refer  to  £c,  which  is  computed 
for  a  background  medium  £m. 

For  a  cylinder  with  an  electrically  small  radius  a,  i.e.,  /3da  «  1 ,  an  approximation 

to  c„  can  be  obtained  by  utilizing  the  small  argument  forms  of  the  Bessel  and  Hankel 
functions  as  follows  [24]: 


(-I)'" 

WH1 

1  fx 

n\l 


1 

n  ! 

UJ 

for  n  <  0 
for  n  =  0 
for  n  >  0 


(3.6) 
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for  n  >  0 


(3.7) 
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If  (J3m  a)  is  also  small,  i.e.  j3ma  «  1 ,  then  the  term  with  index  n  =  0  dominates  the 
series  in  (3.1).  Substituting  (3.6)-(3.9)  into  (3.2),  Co  can  be  approximated  by 

-n(  pma^(£c-\))2 


4  i 


(3.10) 


and  E,  becomes 


4i 


(3.11) 


Figure  3.2  compares  scattered  field  intensities  computed  via  (3.1)  and  (3.11)  for 
cylinders  with  different  radii  and  relative  permittivities  (£c).  As  expected,  the 
approximate  solution  works  well  when  j3da  «  1  and  j3ma  «  1.  In  Figure  3.2,  we  show 
the  results  for  relative  permittivity  ec  <1.5 .  It  is  found  that  the  approximate  solution 
agrees  with  the  exact  solution  for  /3„,  a  <0.15. 
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(a) 


0  0.05  0.1  0.15  0.2  0.25  0.3  0.35 


(b) 

Figure  3.2.  Comparisons  between  scattered  field  intensity  computed  by  the  exact 
solution  and  the  small  radius  approximation  at  p  =  X  and  (j)  =  180°.  Vertical  arrows 
indicate  the  point  at  which  fa  a  =  0.1.  TMZ  case,  (a)  Magnitude,  (b)  Phase. 
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From  (3. 11),  we  can  write 


E,\=xfc\HS>(fi„p)\[aJ\e^i 


(3.12) 


We  find  that  the  result  is  proportional  to  the  cross-sectional  area  of  the  cylinder 
(expressed  in  square  wavelengths)  and  to  the  dielectric  contrast.  If  we  define  an  effective 

radius  ae,TM  =  a->J~  £c  -  1  and  let  J3X  =  max(  J3W ,  Pd ),  then  all  the  results  in  Figure  3.2  can  be 


summarized  by  a  single  universal  curve  shown  in  Figure  3.3.  If  the  medium  permittivity, 
operating  frequency,  radius  and  cylinder  permittivity  are  known,  the  scattered  fields  can 
be  found  from  Figure  3.3  directly.  Note  that  the  universal  curve  is  valid  when  pxa  «  1. 


Figure  3.3.  Scattered  field  intensity  at  distance  p  =  X  as  a  function  of  effective  radius  for 
TMZ  polarization. 

Furthermore,  the  scattered  field  intensity  can  also  be  plotted  as  a  function  of  £c  as 
shown  in  Figure  3.4.  This  figure  shows  that  if  the  cylinder’s  relative  permittivity 
approaches  unity,  the  scattered  fields  will  vanish,  as  expected. 
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Figure  3.4.  Scattered  field  intensity  at  distance  p  =  A  as  a  function  of  £c  for  TMZ 
polarization. 


3.2  Scattering  from  a  Dielectric  Circular  Cylinder:  TE-  Polarization 

The  foregoing  analysis  is  now  repeated  for  a  TEZ  plane  wave  normally  incident  on 
a  dielectric  cylinder  with  radius  a ,  permittivity  £d  and  permeability  p0  (see  Figure  3.5). 


Figure  3.5.  TEZ  plane  wave  impinges  on  a  dielectric  circular  cylinder. 
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For  this  geometry  the  scattered  fields  is  [23] 


where 


Es=pEp  +  0E , 
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Eo  =  magnitude  of  incident  fields  (set  to  1). 
1Jm 


(3.13) 


(3.14) 

(3.15) 

(3.16) 


For  the  electrically  small  cylinder  case  (i.e.,  /3da  «  1  and  jBma  «  1),  the  small 

argument  approximations  of  the  Bessel  and  Hankel  functions  given  in  (3.6)-(3.9)  permits 
us  to  simply  the  coefficients/,,  as 
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(3.17) 


Thus,  for  a  small  radius  cylinder,  Ep«  £/,  and  the  first  two  terms  (n=  +1  and  -1) 
dominate  the  series  of  (3.14),  so  E  becomes 


(~i)7T 


(J3mp) 


P,na, 


£  +1 


(3.18) 
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Figure  3.6  compares  scattered  field  intensities  computed  via  (3.13)  and  (3.18)  for 
cylinders  with  different  radii  and  £c.  As  expected,  the  approximate  solution  works  well 
when  /3d  a  «1  and  J3W  a  «1.  Figure  3.6  shows  that  for  TE  case  with  relative 
permittivities  £c  <  7.5 ,  the  approximate  solution  agrees  with  the  exact  for  [3m  a  <  0.3. 
From  (3.18)  one  can  write 


—  i  xfi2 
E  = — — ^ 


(fi.P) 


n  2 


£  -1 


£  +1 


(3.19) 


Once  again,  we  can  express  these  results  using  a  single  universal  curve  (Figure  3.7)  if  the 
effective  radius  is  defined  as  aejE  =a  S  £c  -1  /  £c  + 1  and  we  let  [3X  =  max(  J3m  ,  Pd  )■ 


Note  that  the  universal  curve  is  valid  when  pxa  «  1. 

Finally,  the  scattered  field  intensity  can  also  be  plotted  as  a  function  of  £c  as 
shown  in  Figure  3.8.  As  expected,  if  the  relative  permittivity  of  the  cylinder  approaches 
unity,  the  scattered  fields  vanish. 
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(b) 

Figure  3.6.  Comparisons  between  scattered  field  intensity  computed  by  the  exact 
solution  and  the  small  radius  approximation  at  p  =  A  and  (j)  =  180°.  Vertical  arrows 
indicate  the  point  at  which  /3da  =  0.1.  TEZ  case,  (a)  Magnitude,  (b)  Phase. 
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Magnitude  of  field  intensity  (dB) 


Figure  3.7.  Scattered  field  intensity  at  distance  p  =  A  as  a  function  of  effective  radius  for 
TEZ  polarization. 


Figure  3.8.  Scattered  field  intensity  at  distance  p  =  / 1  as  a  function  of  £c  for  TEZ 
polarization. 


3.3  Comparisons  of  Magnitude  of  Scattered  Fields  from  TMZ  and  TE-  Cases 

We  now  use  the  foregoing  approximations  to  predict  GPR  backscatterer  from  a 
long,  thin  linear  root.  The  GPR  is  assumed  to  be  far  enough  from  the  root  to  approximate 
a  locally  planar  incident  field.  The  root  radius  is  assumed  to  satisfy  the  small  radius 
approximation,  fy  a  «1,  and  is  modeled  by  an  infinitely  long  cylinder.  Under  these 
conditions,  (3.11)  and  (3.17)  may  be  used. 

Figure  3.9  shows  the  relative  magnitude  of  the  scattered  fields  for  parallel  and 
perpendicular  polarizations  as  a  function  of  relative  permittivity  with  respect  to  the 
medium  (£c).  It  can  be  seen  from  these  figures  that  as  £  increases,  fields  backscattered 
by  the  perpendicular  polarization  can  be  neglected,  which  means  that  only  roots  in  the 
same  direction  as  the  incident  field  need  to  be  considered.  When  £c  is  relatively  low, 
scattering  from  roots  is  less  sensitive  to  both  the  incident  polarization  and  the  root  radius. 
The  intensity  of  the  backscattered  field,  however,  is  relatively,  small  in  that  case. 


Figure  3.9.  Relative  magnitude  of  scattered  fields  for  the  TMZ  and  TEZ  polarizations  as  a 
function  of  relative  permittivity  and  radius. 
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3.4  Summary 

This  chapter  has  reviewed  the  exact  solutions  for  scattering  from  an  infinitely 
long  dielectric  circular  cylinder.  The  small-radius  approximations  for  plane-wave 
scattering  from  a  cylinder  for  both  TMZ  and  TEZ  polarizations  were  described.  The 
approximate  solutions  show  how  the  scattered  fields  depend  on  the  cylinder’s  radius  and 
permittivity.  It  is  also  shown  that  for  a  root  with  a  low  relative  permittivity  the 
backscattered  field  is  less  dependent  on  the  root  orientation,  although  backscattering  is 
much  weaker  in  that  case. 
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4.0  Implementation  of  DDA  to  Model  the  GPR  System 


In  Chapter  2  it  was  shown  that  the  DDA  code  can  accurately  calculate  scattering 
from  cylinders.  In  this  chapter  we  use  the  DDA  to  simulate  GPR  data.  Experiments  and 
data-processing  procedures  for  model  validation  are  described.  Comparisons  of  measured 
and  calculated  data  are  also  presented. 


4.1  Modeling  of  the  GPR  System 

The  GPR  system  used  in  the  experiments  was  described  previously  in  Part  1  of 
this  report  [1],  and  is  illustrated  in  Figure  4.1(a)  and  4.1(b).  The  radar  functions  are 
performed  by  a  network  analyzer  controlled  by  a  laptop  computer.  The  system 
components  are  mounted  on  a  low,  wheeled,  plastic  cart.  The  antenna  is  a  broadband, 
dual-polarization,  dielectric-loaded  horn  antenna  developed  at  OSU/ESL  [25],  and  is 
further  described  below.  The  antenna  and  its  dipole  elements  are  shown  in  Figure  4.1(c). 


(a)  (b) 


Conductor 


Figure  4.1.  (a)  GPR  system  from  top  view,  (b)  GPR  system  from  side  view,  (c) 
Dielectric-loaded  horn  antenna 


29 


Reflections  from  the  air-ground  interface  are  typically  strong  and  can  overwhelm 
near-surface  target  returns.  To  reduce  these  surface  reflections,  the  horn  section  of  the 
antenna  is  filled  with  an  artificial  dielectric  having  a  dielectric  constant  of  five.  If  the  soil 
permittivity  is  also  near  five,  the  antenna  is  matched  to  the  ground,  and  the  reflection 
coefficient  is  minimized.  This  fact  simplifies  modeling  of  the  GPR,  since  the  air-ground 
interface  can  be  removed  as  shown  in  Figure  4.2. 

A  simple  model  for  the  GPR  antenna  is  used  here.  The  antenna  is  represented  by 
two  perpendicular  dipole  elements  oriented  in  the  x  and  y  directions  (assuming  that  the 
interface  lies  in  the  x-y  plane).  Note  that  the  pattern  of  the  real  antenna  may  be  different 
from  a  dipole,  but  it  will  be  shown  later  that  this  simplification  still  works  quite  well.  A 
complete  model  of  this  antenna  on  an  air-ground  interface  using  three-dimensional  finite- 
difference  time  domain  (3-D  FDTD)  has  been  presented  by  Kwan  et  al.  [26]. 


(a)  (b) 

Figure  4.2.  (a)  Geometry  of  the  GPR  system  over  the  soil,  (b)  GPR  system  when  the  air- 
ground  interface  is  removed. 


30 


4.2  Verification  of  the  Model 

Experiments  were  performed  to  verify  the  above  model.  We  first  examined  the 
response  as  a  function  of  frequency  and  distance  from  a  dielectric  rod  of  length 
five  feet  and  diameter  one  inch  with  a  dielectric  constant  of  30.  (The  “rod”  is  actually  a 
hose  filled  with  the  liquid,  methanol.)  The  rod  was  buried  at  a  depth  of  one  foot  in  a  pit 
filled  with  sand  located  at  ESL.  The  GPR  was  scanned  over  the  rod  along  a  path 
perpendicular  to  the  rod  axis.  The  antenna  was  oriented  such  that  the  two  dipole  elements 
were  parallel  and  perpendicular  to  the  direction  of  the  scan.  Those  elements  are  referred 
to  as  Channel  1  and  Channel  2,  respectively.  The  data  were  collected  for  frequencies 
from  0.05  to  1.5  GHz  at  4.875  MHz  increments.  Position-scan  data  were  acquired  at  3 
inch  increments  over  a  total  range  of  72  inches  as  shown  in  Figure  4.3.  A  comparison  of 
measured  and  modeled  time-domain  responses  for  Channel  2  (parallel  to  the  rod)  is 
shown  in  Figure  4.4. 
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Figure  4.3.  The  initial  experimental  set-up.  (a)  Top  view,  (b)  Side  view 
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(b) 

Figure  4.4.  Comparison  of  the  time  domain  responses  of  a  1”  diameter  rod. 
(a)  Measured  data,  (b)  DDA  calculated  data. 
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The  measured  and  calculated  responses  have  similar  temporal  and  spatial 
characteristics  (cf.  Figures  4.4(a)  and  4.4(b),  respectively),  but  the  former  one  has  a  lower 
signal-to-clutter  ratio.  This  clutter  is  a  result  of  soil  inhomogeneities  produced  by  digging 
and  burial  activities. 

In  order  to  obtain  responses  with  reduced  clutter,  a  second  approach  was  used  in 
subsequent  measurements.  The  new  experiment  setup,  shown  in  Figures  4.5  and  4.6, 
permits  us  to  remove  the  majority  of  the  clutter  via  background  subtraction.  A  hole  was 
made  in  a  side  of  the  test  pit  wall  such  that  a  rod  could  be  inserted  from  the  side  without 
disturbing  other  parts  of  the  soil.  After  acquiring  the  response  of  the  background  soil  in 
the  absence  of  the  rod  (see  Figure  4.6(a)),  rods  with  different  dielectric  constants  and 
diameters  were  inserted  into  the  hole  for  measurement  one  at  a  time.  The  data  were 
collected  from  the  polarimetric  GPR  antenna  located  directly  above  the  rod.  Figure  4.6(b) 
shows  a  diagram  of  the  experiment  after  a  rod  was  inserted. 


Figure  4.5.  The  set-up  of  the  experiment. 
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Wall 


(a) 


Figure  4.6.  The  set-up  of  the  experiment,  (a)  GPR  collects  data  from  a  background,  (b) 
GPR  collects  data  from  a  background  and  a  cylinder. 


The  background  subtraction  process  is  straightforward.  Define  the  frequency 
responses  and  transfer  functions  as  follows: 

Hp(0))  =  Transfer  function  of  the  antenna  when  it  radiates. 

B(  co)  =  Frequency  response  of  the  background  soil  (including  the  interface) 
collected  before  the  rod  is  inserted. 

Hb(co)  =  Transfer  function  of  the  antenna  when  it  receives.  (For  a  reciprocal 
antenna,  HB(co)  =  H/i  co). ) 

T(  co)  =  Frequency  response  of  fields  backscattered  from  the  rod. 

CO  =  Angular  frequency. 

In  the  experiment  the  frequency-domain  data  were  “S”  parameters  whose 
definition  is  [V  Output ]  /  [  V  Input ]  [27].  With  the  above  definition,  we  can  express  the 
measured  responses  of  the  background  (  Rj(co)  )  and  background  with  a  rod  (  /Gi CO)  ). 
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R  m  =  [0l,tpun]  =  H B ( (d)Hf  (a)  B{co) 


(4.5) 


R  m  =  VOulpum  =  H B  (co)  H F  (co)  (  B(co)  +  T(co) )  (4.6) 

[Input] 

where  we  have  assumed  that  there  is  no  interaction  between  the  background  soil  and  the 
rod.  The  differences  of  the  measurements  in  then 

Rsub(<o)  =  R2(co)-Rl(co)  =HB(co)HF(co)  T(co) .  (4.7) 

Figures  4.7,  4.8  and  4.9  show  examples  of  the  collected  and  subtracted  data  for 
Ri(co),  R2( (O)  and  Rsub( OJ),  respectively. 
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Figure  4.7.  Frequency  response  of  the  background  in  the  absence  of  a  rod  (R](  (D)). 
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Figure  4.8.  Frequency  response  of  the  background  and  a  rod  ( R2( CO)). 


Figure  4.9.  Difference  between  R2( CO)  and  R/( (d)  (. Rsub( CO)). 


From  the  above  results,  it  is  found  that  R/( 00)  and  R2( OJ)  arc  dominated  by  the 
responses  from  the  antenna  feed  and  the  soil  surface.  As  a  result,  the  subtraction  yields  a 
very  weak  signal  (approximately  -30  dB  below  the  raw  data).  Note  that  the  fast  variation 
in  the  data  shown  in  Figures  4.7  and  4.8  is  caused  by  the  interference  between  the 
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multiple  reflected  signals  in  the  cable  connected  between  the  network  analyzer  and  the 
antenna. 

In  practice,  the  above  subtraction  is  not  perfect  because  there  are  variations  in  the 
system  over  time.  To  reduce  this  effect  Rsub(  0))  is  inverse  Fourier  transformed  to  the  time 
domain,  and  the  residue  from  the  subtraction  at  early  time  (due  to  the  antenna  feed  and 
the  background)  is  gated  out  as  shown  in  Figure  4.10.  The  remaining  signal  is  then 
Fourier  transformed  back  to  the  frequency  domain  to  obtain  a  better  estimate  for 
H b(q))H F((0)  T(co) . 

From  (4.7)  it  follows  that  if  H B{co) H F(co)  is  known,  T(cd)  can  be  found,  and  a 
direct  comparison  with  the  DDA  code  is  possible. 


x  104 


Figure  4.10.  Subtracted  data  in  the  time  domain  and  a  late-time  window  to  remove  early 
time  clutter 

4.3  Finding  the  Antenna  Response  H B  ((0)  H F  (CO) 

An  additional  experiment  must  be  done  to  obtain  the  antenna  response.  The  setup 
and  procedure  are  similar  to  those  shown  in  Figures  4.6(a)  and  4.6(b),  except  that  a  long, 
thin  wire  was  used  as  the  target.  The  frequency  response  of  the  wire  is  known  to  be 
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approximately  independent  of  frequency  (shown  in  Appendix  A),  which  permits  easy 
estimation  of  the  antenna  response. 

The  background  measurement  was  done  exactly  the  same  as  shown  in  Figure 
4.6(a).  For  the  response  of  the  background  and  thin  wire,  the  latter  was  placed  between 
the  antenna  and  the  ground  so  that  the  angle  between  the  thin  wire  and  each  of  the 
antenna  elements  is  45  degrees. 

If  coupling  between  the  ground  and  wire  is  ignored,  subtraction  yields 

Rsub,w(G))  =  HB(0))HF(G))W{0) )  (4.8) 

where  W( CO)  (equal  to  a  constant  C)  is  the  frequency  response  of  a  thin  wire.  When 
absolute  measurements  are  not  being  performed,  C  can  be  set  to  be  1.  The  result  is  shown 
in  Figure  4.11.  By  dividing  (4.7)  by  (4.8),  the  frequency  response  of  the  rod  T( CD)  can  be 
obtained  as  shown  in  Figure  4.12. 


Figure  4.11.  Response  from  a  long  thin  wire  at  the  surface,  which  is  proportional  to  the 
antenna  response  HB(  CO)HF(  CO). 
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Figure  4.12.  The  frequency  response  of  the  rod  T(  CO),  obtained  by  applying  the  window 
in  the  time  domain  and  dividing  by  the  antenna  response. 

4.4  Comparison  of  Experiment  and  DDA  Predictions  for  T(  CO) 

Figure  4.13  compares  the  DDA-calculated  rod  responses  with  T(  to)  acquired  via 
the  above  procedure  for  two  rods  of  the  same  dielectric  constant  (30.0),  but  different 
diameters  (0.5  and  1.0  inches).  Figure  4.14  shows  a  similar  comparison  for  two  rods  with 
the  same  diameter  (1  inch),  but  different  dielectric  constants  (2.4  and  30.0).  In  these 
results  the  calculated  data  have  been  weighted  by  a  constant  factor  that  accounts  for 
antenna  gain  and  RF  losses.  (The  same  factor  is  used  in  all  of  these  results.)  The  rods 
with  dielectric  constants  of  2.4  and  30.0  are  plastic  and  methanol,  respectively.  All  of  the 
rods  are  shown  in  Figure  4.15.  The  above  comparison  shows  good  agreement  between 
the  model  and  the  measurement  although  there  are  some  discrepancies  which  arise 
because  of  three  factors.  First,  the  window  applied  in  the  time  domain  to  remove  the 
early-time  clutter  causes  some  distortion  at  both  ends  of  T(  CO)  as  can  be  seen  in  Figure 
4.13(a)  and  4.14(a).  Second,  the  dielectric  constant  of  the  hose  we  used  differs  from  that 
of  methanol.  This  can  induce  errors  in  the  magnitude  of  T( CO).  Finally,  small  timeshifts  in 
Figure  4.13(b)  and  4.14(b)  between  the  measured  and  calculated  data  may  arise  because 
the  dielectric  constant  of  the  soil  may  not  be  exactly  5. 
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Figure  4.13.  Comparison  of  measured  and  DDA-calculated  backscatterer  for  two  rods 
with  the  same  permittivity  and  different  diameters,  (a)  Frequency  response,  (b)  Time- 
domain  response. 
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(a) 


Time  (ns) 


(b) 

Figure  4.14.  Comparison  of  measured  and  DDA-calculated  backscatterer  for  two  rods 
with  the  same  diameter  and  different  permittivities,  (a)  Frequency  response,  (b)  Time 
domain  response. 
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Figure  4.15.  Dielectric  rods  used  in  the  experiment.  Note  that  “d”  stands  for  diameter. 
4.5  Summary 

The  measurement  approaches  used  at  the  OSU/ESL  test  site  for  collecting  wide 
band,  dual-polarized  data  of  buried  rods  made  of  plastic  and  methanol  have  been 
discussed.  Both  position- scanned  and  single-position  data  were  collected.  Burying  the 
rods  disturbed  the  soil,  which  produced  clutter  of  a  magnitude  similar  to  the  dielectric  rod 
return.  Single-position  measurements  were  found  to  provide  better  data  for  model 
validation.  Comparisons  of  calculated  and  measured  returns  from  buried  rods  show  good 
agreement.  This  demonstrates  that  the  DDA  code  can  be  used  to  accurately  compute  the 
fields  scattered  from  a  root  structure  for  a  wide  range  of  dielectric  constants  and 
diameters. 
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5.0  Case  Studies  of  Effects  of  Roots  on  GPR 


The  simple  model  developed  in  Chapter  4  for  simulating  GPR  scans  will  be  used 
to  study  scattering  from  the  root  structure  shown  in  Figure  2.1(b).  This  chapter  begins 
with  a  discussion  of  how  root  scattering  varies  as  a  function  of  root  dielectric  constant 
and  diameter.  Comparisons  between  experimental  and  simulated  GPR  scans  will  be 
presented.  Two  issues  concerning  the  effects  of  roots  on  GPR  system  are  studied,  namely, 
the  increase  in  clutter  level  and  root  attenuation  of  targets  signatures. 


5.1  Properties  of  Root  Scattering 

Before  presenting  a  comparison  of  measured  and  modeled  results,  we  first  use  the 
model  to  illustrate  the  dependence  of  root  scattering  on  dielectric  constant  and  frequency. 
A  long,  thin,  isolated  root  will  be  examined.  The  results  can  be  contrasted  with  those  in 
Chapter  3  for  the  thin,  infinitely  long  root.  The  configuration  for  this  case  is  shown  in 
Figure  5.1(a)  and  5.1(b).  Note  again  that  the  radiating  source  is  modeled  as  two 
perpendicular  short  dipole  elements  oriented  in  the  x  and  y  directions  (assuming  that  the 
interface  lies  in  the  x-y  plane).  Throughout  this  chapter,  we  will  assume  the  dielectric 
constant  of  roots  (£ct)  to  be  20+10i  which  is  calculated  using  the  Debye  Cole  dual¬ 
dispersion  model  described  in  Appendix  B.  The  soil  dielectric  constant  and  conductivity 
are  5  and  0.0035  S/m,  respectively,  as  measured  by  the  Army  Research  Laboratory 
(ARL).  The  details  are  explained  in  Appendix  C. 


GPR  Antenna 
Soil 


root 


(a) 


Soil,  £jn  —  5, 

<7=  0.0035  S/m 


Center  of  root 


1.5  m- 

(b) 


Figure  5.1.  Configuration  of  the  GPR  scan  used  to  study  the  scattering  characteristics  of 
a  long,  thin  root  with  dielectric  constant  20+10i.  (a)  Front  view,  (b)  Side  view. 
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The  fields  scattered  from  the  root  are  computed  at  the  antenna  and  are  normalized 
by  the  incident  field  that  would  exist  in  the  absence  of  the  root  at  the  location  of  the 
root’s  centroid.  Figure  5.2  shows  the  calculated,  scattered  electric  field  intensity  for  the  x 
polarization  as  a  function  of  both  diameter  measured  in  wavelengths  (D/Ad)  and  root 
dielectric  constant. 


005  o.i  tr.is  02  o.;b 


D  /  Ad 

Figure  5.2.  Magnitude  of  the  scattered  field  intensity  (in  dB)  as  a  function  of  root 
diameter  in  wavelengths  and  dielectric  constant. 

As  expected,  the  magnitude  of  the  field  scattered  from  a  long  thin  root  increases 
as  the  diameter  increases.  It  also  increases  with  the  dielectric  contrast  £  =  £c\/£m  where  £m 
is  the  dielectric  constant  of  the  soil.  These  simple  dependencies  on  £c  - 1  and  D/Ad 
were  also  found  in  Chapter  3  for  an  infinitely  long  cylinder  described  by  equation  (3.11). 
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5.2  Approximations  to  Reduce  Model  Complexity 

It  was  observed  during  early  tests  of  the  DDA  code  that  taking  into  account  all  the 
roots  in  Figure  2.1(b)  at  every  antenna  position  is  computationally  expensive,  because  of 
the  large  number  of  unknowns.  Fortunately,  further  approximations  can  be  employed  to 
reduce  the  size  of  the  problem. 

In  general,  soil  is  lossy  and  the  antenna  has  mildly  directional  radiating  and 
receiving  patterns.  This  implies  that  roots  far  from  the  antenna  and  outside  the  main 
pattern  will  have  little  contribution  to  the  scattered  fields.  Thus,  instead  of  considering  all 
roots  at  each  antenna  position,  it  is  sufficient  to  consider  only  those  roots  that  both  lie 
within  some  distance  of  the  source  and  are  within  the  transmitting/receiving  patterns  of 
the  antenna.  In  this  work  we  specify  that  the  soil  has  a  dielectric  constant  of  5, 
conductivity  of  0.0035  S/m,  and  that  the  antennas  arc  described  by  dipoles.  As  a  result, 
only  scattered  fields  caused  by  root  segments  within  the  cone  shown  in  Figure  5.3  need  to 
be  considered  in  the  calculation.  The  angle#  is  found  to  be  76.8°,  which  is  defined  by  the 
criterion  that  the  incident  field  intensity  at  Point  P  is  at  least  25  dB  below  that  at  Point  Q. 


Ground 

Surface 


Figure  5.3.  Only  roots  in  a  circular  cone  shape  are  considered  in  the  scattering  compu¬ 
tation. 


Another  simplification  involves  the  condition  in  equation  (2.25),  namely, 


^(2«max)  - 1- 


(5.1) 


The  above  condition  leads  to  an  overly  fine  sampling  of  the  problem  when  a  large 
number  of  roots  arc  taken  into  account,  and  it  is  reduced  to 


k  ^max  ^  1  • 


(5.2) 
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Figure  5.4  shows  good  agreement  between  computed  scattered  field  intensities  as 
a  function  of  frequency  after  the  above  approximations  are  applied.  Here,  the  depth  of  the 
roots  is  6  inches,  and  their  dielectric  constant  is  20+10i.  The  antenna  is  positioned  above 
the  large  root  mass  in  Figure  2.1(b).  This  approach  reduced  the  number  of  unknowns  at  1 
GHz  from  31,152  to  9,960,  which  makes  the  calculation  more  efficient. 


Figure  5.4.  Comparison  of  scattered  field  intensities  computed  when  considering  all 
roots  discretized  using  Equation  (5.1)  and  only  those  roots  that  fall  within  the  circular 
cone  discretized  using  Equation  (5.2). 
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5.3  Comparison  Between  Modeled  and  Measured  GPR  Data 

To  demonstrate  that  the  model  can  reproduce  fields  scattered  from  a  realistic  root 
system,  model  results  were  compared  with  data  collected  at  Eglin  AFB  along  the  path 
shown  in  Figure  5.5.  Since  it  was  shown  in  the  first  report  [1]  that  the  radar  is  insensitive 
to  the  roots  for  frequencies  higher  than  1.0  GHz,  we  will  only  consider  the  frequency 
range  from  0.05  to  1.0  GHz.  The  frequency  step  for  the  measured  data  is  4.875  MHz  with 
three  inch  increments  in  position,  but  for  the  simulated  data,  the  frequency  step  is  50 
MHz  with  0.1333  m  increments.  Note  that  in  the  simulation,  the  depth  of  the  root 
structure  is  set  to  be  6  inches.  Figure  5.6  shows  the  calculated  and  measured  data  plotted 
as  functions  of  time  delay  and  scan  position.  Within  the  modeled  interval,  similarities 
between  the  measurement  and  model  can  be  seen.  The  arrows  indicate  to  the  positions  at 
which  strong  responses  are  expected  and  they  agree  well  with  the  root  location  shown  in 
both  figures.  Some  disagreement  between  the  measurement  and  model  can  be  attributed 
to  our  two-dimentional  (2-D)  root  model:  the  true  root  structure  has  some  variation  in 
depth,  but  the  model  places  all  roots  at  a  fixed  depth. 
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Figure  5.5.  GPR  scan  path  and  root  structure 
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Figure  5.6.  Comparison  of  measurements  and  model  predictions,  (a)  Experimental  data, 
(b)  Model  simulation. 
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5.4  Case  Studies  of  the  Effects  of  Roots  on  GPR  Systems 

From  the  above  comparison,  it  is  evident  that  the  model  can  predict  root  scattering 
fairly  well.  Hence,  the  model  can  be  used  with  some  confidence  to  study  the  effects  of 
roots  on  GPR  performance.  We  begin  with  a  study  of  two  issues:  root  generated  clutter 
and  the  tendency  of  roots  to  produce  target  attenuation. 


5.5  Target  Detection  in  the  Presence  of  Roots 

We  now  evaluate  the  detectability  of  a  known  UXO  target  in  a  clutter  field  caused 
by  root  scattering.  A  complete  analysis  of  this  topic  is  beyond  the  scope  of  this 
investigation,  but  useful  insight  can  be  drawn  from  a  relatively  simple  analysis. 

It  is  assumed  that  the  radar  emits  a  (complex)  waveform /t(Y).  We  further  assume  that 
the  target  response  is  dominated  by  the  specular  return,  which  can  be  expressed  as  Ah(t-Ti),  in 
which  T]  is  a  time  delay  and  A  is  a  complex  scalar  that  accounts  for  both  propagation  losses 
and  a  possible  change  in  phase  due  to  reflection.  The  root  return  is  also  approximated  by  its 
specular  contribution,  which  we  write  as  Bh(t-T2)  in  which  z?  and  B  are  an  appropriate  phase 
shift  and  a  complex  scalar,  respectively.  The  effects  of  receiver  noise  are  assumed  to  be 
insignificant  in  comparison  to  the  clutter.  The  detection  problem  can  be  expressed  as 
follows: 


x(t)  = 


j Ah(t  —  T,)  +  Bh(t  -  t2  ), 
|  Bh(t-r2 ), 


if  Hx 

if 


(5.3) 


in  which  x(t)  is  the  measured  signal,  and  Ho  (H f)  is  the  hypothesis  that  the  target  is  absent 
(present)  in  the  measurement. 

An  appropriate  technique  for  detecting  such  signals  involves  the  matched  filter,  in 
which  the  received  waveform  is  convolved  with  the  time-reversed,  conjugated, 
transmitted  waveform.  We  find 


x(t)®h\-t) 


J Av(t  -  Tj )  +  Bv(t  -  t2),  if  H, 
{  Bv(t  —  T2 ),  if  H0 


in  which  v(t)  is  the  response  of  the  matched  filter  to  the  transmitted  waveform 

v(t)  =  h(t)  ®h*(—t)  . 


(5.4) 


(5.5) 
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A  practical  implementation  of  this  concept  would  involve  searching  the  matched 
filter  output  for  local  maxima  of  the  magnitude  (which,  presumably,  would  be  the 
maxima  of  v(t)).  A  detection  decision  would  be  required  for  each  maxima. 

If  the  roots  and  targets  are  always  found  at  different  depths,  then  detection  can  be 
as  simple  as  range-gating  the  matched  filter  output  and  applying  a  threshold.  The  case  of 
interest  here,  however,  is  when  both  the  target  and  roots  are  in  the  same  range  bin.  In 
that  case,  we  can  take  Ti~T2.  If  Vo  is  the  maximum  value  of  v(t),  then  the  detection 
problem  is  equivalent  to  a  test  of  the  scalar 

max[|  x(t)  ®h*  (—t)  |]  f|  A  +  B 
v0  \  |fl|, 
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Figure  5.7.  Sixteen  chosen  root  locations. 
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Figure  5.8.  Sixteen  examples  of  clutter  spectra  computed  over  known  root  locations. 


When  the  probability  densities  of  A  and  B  are  known,  then  one  can  formulate  an 
optimal  detection  test  for  this  problem.  In  the  present  case,  however,  the  available 
information  is  limited,  and  less  rigorous  assessment  is  appropriate. 

Consider  first  the  density  of  \B\.  Examples  of  clutter  spectra  for  16  locations 

chosen  directly  over  roots  (see  Figure  5.7)  appear  in  Figure  5.8.  These  data  indicate  that 

the  clutter  is  strongly  non-Gaussian.  On  the  logarithmic  (dB)  scale  used  here,  there  is 
some  degree  of  central  tendency,  and  it  is  plausible  to  approximate  the  probability  density 
values  with  a  log-normal  density  f(u)  given  by 

f(u)du  =  J_  e-<Mu)-M?l(2a2]du  (5.7) 

V  Ixoii 

where  jLl  and  <7 are  constants.  Using  the  change  of  variable 

v  =  ln(«)  (5-8) 


we  obtain  the  alternate  form 
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fiu)du  =  _^g-(v-/b2/(2cF)  dv 

yjl7TO 


(5.9) 


which  shows  that  the  logarithm  of  the  independent  variable  u  has  a  normal  density  in  the 
transformed  variable  v  with  mean  //  and  standard  deviation  <7.  If 

20  v 


5  =  20  log  10(h)  = 


ln(10) 


(5.10) 


is  the  value  of  the  clutter  return  expressed  in  dB,  then 

1 


f(u)du 


4l7 To' 


e^~M'fK2^)ds 


(5.11) 


which  is  a  standard  normal  density  with  the  following  standard  deviation  and  mean  (both 
expressed  in  dB) 

o'  =  20o71n(10) 

//'  =  20///ln(10) 


(5.12) 


These  quantities  can  be  computed  directly  from  the  data  expressed  in  dB. 

Confidence  intervals  for  this  density  are  readily  obtained  from  those  for  a 
standard  normal  density,  and  they  permit  us  to  define  appropriate  detection  thresholds  for 
a  Neyman- Pearson  receiver  [28]  (i.e.,  one  which  achieves  a  specified  probability  of  false 
alarm.)  For  this  receiver,  knowledge  of  the  density  of  A  is  not  required  to  estimate  the 
probability  of  false  alarm.  (They  are,  however,  necessary  to  estimate  the  probability  of 
detection.)  For  example,  if  the  detection  threshold  for  the  target  density  in  dB  is  set  at 
ju'+1.28(f,  then  the  probability  of  false  alarm  is  0.1.  At  a  threshold  of  ju'+2.33&  the 
false  alarm  probability  is  reduced  to  0.01.  The  relevant  false  alarm  rate  is  obtained  by 
dividing  these  probabilities  by  the  sensor’s  pixel  area. 

Use  of  the  foregoing  results  is  illustrated  in  Figure  5.9(a),  in  which  we  show  a  105 
mm  shell  signature  at  depth  0.35  m  computed  for  several  shell  angles  via  FDTD  [26]. 
Figure  5.9(b)  defines  the  shell  angle.1  This  target  signature  is  superimposed  with  the 
average  clutter  spectra  of  Figure  5.8.  The  decision  threshold  (a  function  of  frequency)  for 
a  false  alarm  rate  of  0.01/m“  (assuming  1  m"  sensor  pixels)  is  also  shown. 


1  The  shell  orientation  in  these  results  is  nose-up,  while  real  UXO  signatures  will  typically  be  found  nose- 
down.  For  the  purposes  of  this  analysis,  however,  the  difference  in  the  target  signature  will  be  small. 
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Figure  5.9.  (a)  Detectability  calculations  for  a  105  mm  shell  in  roots.  The  indicated  line 
is  the  decision  threshold  for  a  false  alarm  rate  of  0.01/m2.  (b)  Orientation  of  105  mm 
shell  as  a  function  of  0. 


The  comparison  in  Figure  5.9(a)  suggests  that  responses  from  the  shell  exceeds  the 
mean  root  clutter  level  for  most  frequencies  and  shell  angles.  Nonetheless,  it  is  only  in  the 
frequency  regime  below  0.5  GHz  that  the  target  signature  exceeds  the  threshold  required 
for  a  10 ’“/nr  false  alarm  rate.  Based  on  these  data,  one  expects  better  performance  for  low- 
frequency  waveforms,  and  significant  challenges  in  achieving  very  low  false-alarm  rates. 
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5.6  Target  Signature  Attenuation 

Next,  consider  the  UXO-shaped  lossy  dielectric  target  shown  in  Figure  5.10  (10 
cm  diameter  by  40  cm  length)  with  dielectric  constant  1+20/  placed  under  the  root  mass 
at  depth  0.3  m.  The  target  orientation  is  parallel  to  the  vertical  axis.  The  antenna  position 
is  directly  above  the  root  mass  and  polarized  in  the  same  direction.  To  isolate  the  effects 
of  root-related  attenuation,  frequency  responses  are  computed  for  the  cases  when  the 
target  is  and  is  not  present.  The  result  of  subtracting  these  two  data  sets  is  compared  with 
the  response  from  the  target  when  no  roots  are  present.  The  comparison  in  the  frequency 
and  time  domains  are  shown  in  Figure  5.11(a)  and  5.11(b).  Note  that  in  the  calculation, 
we  include  the  interaction  between  the  roots  and  target.  In  these  results  it  is  important  to 
disregard  the  relative  magnitudes  of  the  target  and  root  responses.  Because  of  limitation 
in  the  DDA  code,  the  target  (a  good  conductor)  was  approximated  by  a  lossy  dielectric 
(£>=1+20/ ),  which  has  an  abnormally  weak  return. 

The  results  suggest  that  roots  have  a  relatively  small  effect  on  the  target  signature. 
They  attenuate  the  target  response  slightly  at  high  frequencies,  and  they  introduce  a  small 
time  delay.  Note  that  the  attenuation  can  reach  3  to  4  dB  at  the  highest  frequency.  These 
effects  can  be  attributed  to  distortion  in  the  incident  field.  Figure  5.12,  5.13  and  5.14  are 
plots  of  the  calculated  field  magnitude  at  different  positions  normalized  by  the  incident 
field  for  a  depth  of  0.3  m.  The  figures  show  the  response  at  the  frequencies  0.1,  0.6  and 
1.0  GHz,  respectively.  The  purple  triangular  shape  illustrates  the  position  of  the  antenna 
which  is  modeled  as  a  dipole  with  polarization  pointing  down  the  page.  It  can  be  seen 
that  at  relatively  high  frequencies,  the  root-induced  distortion  of  the  incident  field  can 
reach  3  dB,  but  for  low  frequencies  this  distortion  is  relatively  small. 


Figure  5.10.  UXO-shaped  target  modeled  with  cylinders 
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Figure  5.12.  Distortion  of  the  incident  field  magnitude  at  the  depth  0.3  m  for  0.1  GHz. 
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Figure  5.14.  Distortion  of  the  incident  field  magnitude  at  the  depth  0.3  m  for  1.0  GHz. 
5.7  Imaging  of  Simulated  GPR  Data 

The  GPR  response  of  a  point  target  extends  for  a  significant  distance  (more 
precisely,  it  produces  a  hyperbolic  arc)  in  the  along-scan  dimension  (more  precisely,  it 
produces  a  hyperbolic  arc),  which  complicates  target  position  estimates.  One  can, 
however,  use  imaging  techniques  to  both  localize  the  target  signature  and  to  increase  the 
signal  to  noise  ratio. 

In  this  section  we  explore  the  effectiveness  of  imaging  in  the  presence  of  root 
clutter.  The  imaging  algorithm  used  here  is  based  on  the  so-called  “frequency- 
wavenumber  ( co-k )  migration”  technique  used  in  exploration  geophysics.  A  detailed 
derivation  of  this  technique  is  explained  in  several  works  [29-31].  For  a  brief  review,  one 
can  refer  to  [32]. 

The  approach  that  we  will  employ  here  is  based  on  the  scalar  wave  equation  and 
the  concept  of  the  “exploding  reflector”  model  [33].  The  latter  can  be  explained  as 
follows.  Assume  that  a  GPR  antenna  collects  data  f(x,z=0,t)  where  x  is  the  scan 
dimension.  The  signal  scattered  back  from  a  target  at  the  antenna  is  the  same  as  that 
received  by  a  passive  sensor  when  (a)  the  target  becomes  a  source  of  radiation  and  (b)  the 
velocity  of  wave  propagation  is  reduced  by  a  factor  of  2. 
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The  spectral  representation  of  the  measured  data  f(x,z=0,t)  can  be  written  as 

F(kx ,  z  =  0,  co)  =  J dx J dt  f{x,  z  =  0,t )  e+i<a~iKx .  (5.13) 

Since  we  use  the  exploding  reflector  model,  the  field  distribution  of  interest  is  f(x,z,t=0). 
It  is  possible  to  show  that  this  distribution  given  by 

f(x,z,t  =  0)  =  —^\dx\dkx  e+ik*xF(kx,z  =  0,6))e+iK*z  (5.14) 

(InY J  J 

where  K_  =  ^jk2  -  k  v  2  and  k=2co/  v.  The  field  distribution  in  Equation  (5.14)  will  be 

plotted  as  an  image.  In  evaluating  the  transforms  in  Equation  5.13,  we  apply  a  weighting 
function,  i.e.,  a  two-dimensional  window,  to  f(x,z=0,t),  which  restrict  the  images  to  the 
center  of  the  field  of  view. 

The  effectiveness  of  imaging  is  illustrated  in  Figure  5.15(a),  which  shows  the 
normalized  time-domain  response  of  a  UXO-shaped  target  buried  at  a  depth  of  0.55  m 
under  the  root  mass.  Figure  5.15(b)  shows  the  normalized  radar  image.  The  bright  spots 
in  the  image  correspond  to  the  positions  of  the  root  mass  and  the  target,  and  it  is  evident 
that  roots  can  produce  strong  image  artifacts.  This  suggests  that  imaging  is  not 
necessarily  an  effective  means  of  clutter  suppression.  Figure  5.16  shows  the  image  of  the 
target  in  the  absence  of  roots.  The  data  were  normalized  by  the  same  factor  used  in 
Figure  5.15(b).  By  comparing  the  above  figures,  it  is  found  that  roots  can  attenuate  the 
relative  magnitude  of  the  image,  but  the  effect  is  quite  small. 

From  the  discretized  root  representation  in  Figure  2.1(b)  it  is  evident  that,  for 
GPR,  roots  are  more  like  discrete  scatterers  than  random  media.  As  a  result,  image  of  the 
roots  may  be  of  interest  also.  An  example  is  shown  in  Figure  5.16,  which  presents  the 
image  of  Figure  5.6(b).  It  can  be  seen  that  the  locations  of  the  bright  spots  indicate  the 
positions  of  roots.  The  absence  of  a  return  near  the  initial  and  final  ends  of  the  scan  is  the 
result  of  the  window  function. 
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Figure  5.15.  (a)  Time  domain  response  of  a  UXO-shaped  lossy  dielectric  target,  (b) 
Image  of  Figure  5.15(a). 
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Figure  5.16.  Image  of  the  lossy  dielectric  target. 
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Figure  5.17.  Image  of  Figure  5.5(b). 
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5.8  Summary 

The  scattering  characteristics  of  a  long,  thin  root  as  a  function  of  its  diameter  in 
wavelengths  {D/Ad)  and  dielectric  contrast  relative  to  the  medium  ( £,  )  have  been  studied. 
It  was  found  that  the  scattered  field  intensity  is  proportional  to  D/Ad  and  |  ec  —  1 1 .  These 

simple  dependencies  were  also  found  in  Chapter  3  for  an  infinitely  long  cylinder. 
Comparisons  of  measured  and  simulated  GPR  scans  showed  that  the  DDA  model  is 
adequate  to  predict  root  scattering.  Some  discrepancies  were  noted  and  have  been 
attributed  to  the  fact  that  the  modeled  roots  were  all  buried  at  a  depth  of  at  6  inches,  but 
the  true  root  structure  has  some  variations  in  depth.  Further  studies  showed  that  roots 
increase  the  clutter  level,  which  can  cause  problems  in  target  detection.  Based  on  the 
simulated  data,  better  performance  can  be  expected  by  operating  the  system  at  lower 
frequencies.  Roots  can  also  distort  targets’  response,  but  the  effect  is  relatively  small. 
Finally,  experiments  with  imaging  suggest  that  roots  behave  more  like  discrete  scatterers 
than  like  a  continuous  random  medium.  As  a  result,  radar  images  show  distinct  returns 
from  roots,  which  are  not  immediately  distinguished  from  targets. 
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6.0  Conclusions  and  Future  Work 


In  this  report  a  root  structure  obtained  during  experiments  performed  at  Eglin 
AFB,  FL  was  studied  in  detail.  These  roots  were  discretized  and  represented  by  an 
ensemble  of  cylindrical  segments  having  variable  radii.  Scattering  from  the  above 
cylinders  was  computed  via  the  discrete  dipole  approximation  (DDA)  algorithm,  which  is 
an  efficient  technique  for  computing  scattering  from  dielectric  bodies.  As  a  test,  our  DDA 
code  was  applied  to  the  problem  of  plane  wave  scattering  from  a  cylinder.  Comparisons 
of  the  radar  cross  section  computed  by  DDA  and  by  a  body  of  revolution  (BOR)  code 
showed  good  agreement. 

The  exact  solutions  for  both  TM  and  TE  plane-wave  scattering  from  an  infinite 
dielectric  cylinder  were  reviewed,  and  small  radius  approximations  were  derived.  The 
approximate  solutions  have  a  simple  form,  which  illustrates  how  scattered  fields  vary  as  a 
function  of  cylinder  diameter  and  frequency.  The  results  indicate  that  the  magnitude  of 
fields  scattered  from  a  cylinder  increases  as  the  cross-sectional  area  of  the  root  and  the 
permittivity  contrast.  We  also  presented  two  so-called  “universal  curves,”  from  which  the 
magnitude  of  the  scattered  fields  can  be  found  directly  for  any  cylinder  diameter  or 
permittivity. 

Use  of  the  DDA  algorithm  to  model  GPR  scans  was  also  discussed.  In  this  work 
the  GPR  is  assumed  to  use  a  surface-contacting  dielectric-loaded  horn  antenna  matched 
to  the  ground,  and  the  air-ground  interface  was  neglected.  To  model  the  antenna,  two 
perpendicular  dipoles  parallel  to  the  surface  were  used  as  radiating  sources.  The  above 
model  was  verified  by  comparing  calculated  scattered  fields  with  measurements  from 
buried  cylinders  having  different  diameters  and  permittivities.  It  was  found  that  they 
agree  well  in  both  the  frequency  and  time  domains. 

After  validation  of  the  model,  scattering  from  roots  was  explored.  We  began  by 
computing  the  scattering  characteristics  of  a  root  as  a  function  of  frequency  and  diameter. 
The  trend  obtained  matched  that  found  for  scattering  from  an  infinite  dielectric  cylinder. 
Comparisons  of  spatially  scanned  time-domain  responses  obtained  from  both 
measurements  at  Eglin  and  model  predictions  were  also  shown  to  be  in  good  agreement. 
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Two  issues  concerning  the  effects  of  roots  on  GPR  performance  were  studied, 
namely,  root-generated  clutter  and  root-related  target  signature  distortion.  GPR 
responses  directly  over  a  number  of  roots  were  computed  and  fitted  to  a  log-normal 
probability  density.  Analysis  showed  that  this  root-related  clutter  would  produce  high 
false  alarm  rates,  particularly  for  sensors  that  had  significant  frequency  content  above 
about  0.5  GHz.  Simulations  of  target  signatures  with  and  without  overlying  roots  showed 
that  the  roots  produce  little  signature  distortion,  adding  only  a  small  amount  of 
attenuation  and  a  modest  time  delay.  Finally,  attempts  to  image  root  signatures  were 
successful,  which  suggest  that  roots  are  best  treated  as  discrete  clutter  artifacts,  rather 
than  a  continuous  random  phase  screen.  In  summary,  we  concluded  that  the  primary 
effect  of  roots  is  to  add  a  randomly  distributed  set  of  discrete  clutter  sources,  and  that  any 
root-related  loss  in  performance  can  be  partially  offset  by  using  lower  frequencies. 

A  number  of  issues  should  be  explored  in  any  extension  of  this  effort.  First,  the 
current  form  of  the  DDA  code  does  not  handle  targets  that  are  perfect  electric  conductors 
(PECs).  It  is  suggested  that  the  code  be  extended  to  include  those  targets. 

Second,  for  a  GPR  antenna  that  does  not  contact  the  surface  such  as  BoomSAR, 
the  air-ground  interface  must  be  included  in  the  model.  To  deal  with  the  interface 
rigorously,  it  will  be  necessary  to  include  Sommerfeld  integrals,  which  are 
computationally  expensive.  A  fast  algorithm  for  computing  those  integrals  would  be 
essential  for  this  effort. 

Third,  additional  root  samples  should  be  investigated.  This  analysis  was  focused 
on  a  single  root  structure,  which  occupied  an  area  of  only  4  m“.  From  this  small  sample  it 
is  difficult  to  identify  and  characterize  root  anomalies  that  drive  system  performance  in 
an  application  such  as  UXO  detection,  which  is  defined  by  the  occurrence  of  low 
probability  events. 

Finally,  the  universal  curves  should  be  extended  to  approximate  the  fields 
scattered  from  an  average  root  structure  such  that  we  can  approximate  average  root 
returns  without  using  the  DDA  code. 
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Appendix  A  :  Frequency  Response  of  an  Infinitely  Long 
Perfectly  Conducting  Wire  in  the  Frequency  Range  from  0.05- 
1.5  GHz 


In  this  appendix  we  show  that  the  frequency  response  of  an  infinitely  long, 
perfectly  conducting  wire  varies  slowly  for  frequencies  in  the  range  0.05  to  1.5  GHz. 
This  result  justifies  our  decision  to  approximate  wire  scattering  as  a  frequency 
independent  constant,  which  was  done  in  our  calibration  procedure. 

Consider  the  problem  of  plane  wave  scattering  from  an  infinitely  long  cylinder 
with  a  radius  a  =  0.5  mm  in  a  medium  with  dielectric  constant  £r  =  5.  The  problem 
geometry  is  illustrated  in  Figure  Al. 


Figure  Al.  TMZ  plane  wave  impinges  on  a  conductive  circular  cylinder. 


It  can  be  shown  that  if  we  excite  the  problem  with  a  TEZ  uniform  plane  wave,  a 
field  scattered  from  a  conductive  thin  cylinder  is  negligible  [23].  As  a  result,  only  the 
TMZ  case  will  be  considered  here.  With  the  suppression  of  the  factor  e',ca,  the  field 
scattered  from  a  TMZ  plane  wave  can  be  expressed  as  [23] 


E 


s 


n—0 


H^(p^.p)cos{n(f)) 


(A- 1) 


where 


Jn(.)  =  the  Bessel  function  of  the  first  kind  of  order  n  (  n  is  an  integer). 


H"\.)  =  the  Hankel  function  of  the  first  kind  of  order  n. 
Eo  =  magnitude  of  incident  fields  (set  to  1). 
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(A-2) 


P 


if  n  =  0 
if  n  ±  0 


(A-3) 


For  an  electrically  thin  wire  the  quantity  /5a  «  1.  Using  the  small  argument 
approximations  to  the  Bessel  and  Hankel  functions  described  in  Chapter  3,  it  is  clear  that 
the  n= 0  term  dominates  the  series  in  (A-l).  We  can  write 
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Figure  A2  shows  the  magnitude  of  the  backscattered  field  as  a  function  of 
frequency.  For  the  entire  frequency  range  shown,  the  variation  in  magnitude  is  about  4 
dB,  but  for  most  of  the  range,  i.e.,  from  0.3  to  1.5  GHz,  the  variation  is  2  dB,  which  is 
quite  small.  As  a  result,  we  can  approximate  the  frequency  response  of  a  long  conducting 
wire  as  independent  of  frequency  with  a  small  error,  which  arises  primarily  at  low 
frequency. 


Figure  A2.  Magnitude  of  the  field  backscattered  from  a  wire  of  radius  0.5  mm  in  a 
medium  with  relative  permittivity  £r=  5  as  a  function  of  frequency.  Here  p  =  0.3  m  and  (j) 
=  180°. 
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Appendix  B:  Dielectric  Constant  of  Tree  Roots 


In  this  work  we  employed  a  root  dielectric  constant  model  known  as  the  “Debye- 
Cole  dual-dispersion  dielectric  model,”  which  was  developed  by  Ulaby  et  al.  [33].  The 
model’s  authors  asserted  that  it  provides  a  good  estimate  of  the  dielectric  constant  for 
leaves,  stalks  and  trunks,  and  herein  we  have  applied  it  to  roots  as  well.  In  this  model  the 
following  three  components  contribute  to  the  dielectric  constant 

a)  dry  bulk  vegetation  material, 

b)  free  water,  and 

c)  bound  water. 

Only  the  formulae  required  to  compute  the  dielectric  constant  are  shown  here.  Further 
details  are  given  in  [33].  The  expressions  for  the  dielectric  constant  of  vegetation  £v  is 

£v  =  £  r  +  vf  £  f  +  vb£b  (B-l) 


where 


£r  =1.7  -0.1 4 Mg  +  6.16 Mg2 
vfiv  =  Mg  (0.55  Mg  -0.076) 


75  /18  c 
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4.64  Mg 2 
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£h  -  2.9  + 


55 


1  + 


\~i)f 

0.18 


.  0.5 


(B-6) 


(j  =  1.27  (B-7) 

Mg  =  gravimetric  moisture  content,  defined  as 

[  weight  of  water  in  a  root  (g)  ]  /  [  weight  of  a  root  (g)  ] 

/  =  frequency  in  GHz. 


Measurements  at  ESL  of  the  gravimetric  moisture  content  for  five  root  samples 
acquired  at  Eglin  yielded  the  Mg  values  0.450,  0.500,  0.568,  0.660,  and  0.730.  The 
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average  value  is  0.58.  For  this  value,  the  real  and  imaginary  parts  of  £,  are  computed  from 
(B-l)  and  shown  in  Figure  B1  as  a  function  of  frequency.  Note  that  this  model  is  valid 
(produces  error  less  than  ±  20%  )  for  the  frequency  range  from  0.5  to  20.0  GHz.  Since  we 
are  interested  in  the  frequency  range  from  0.05  to  1.0  GHz,  we  will  not  use  the  exact 
model  in  the  DDA  code,  but  the  round  number  of  £v  instead.  Throughout  the  report,  £v  is 
set  to  be  20+10z. 


Figure  Bl.  Real  and  imaginary  parts  of  £,  of  a  root  with  gravimetric  moisture  content  = 
0.58. 
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Appendix  C:  Dielectric  Constant  of  Soils 

Samples  of  soils  from  Eglin  AFB,  FL  were  collected  at  two  sites  by  a  team  from 
ESF.  As  described  in  the  first  report  [1],  these  sites  were  chosen  to  have  high  and  low 
root  densities,  and  are  referred  to  as  “Site  1”  and  “Site  2,”  respectively.  Soils  were 
collected  at  six  inches  depth  intervals  from  surface  to  three  feet.  The  dielectric  constant 
and  conductivity  of  samples  from  Site  1  were  measured  by  the  U.S.  Army  Research 
Faboratory  (ARE). 

Soil  moisture  contents  from  Site  1  and  Site  2  were  measured  at  ESF,  and  the  data 
are  shown  in  Table  Cl  and  C2,  respectively.  The  data  are  plotted  in  Figure  Cl.  Figures 
C2,  C3  and  C4  show  the  dielectric  constant  and  conductivity  of  the  soil  as  a  function  of 
its  depth  and  frequency  (0.1-3  GHz)  for  soil  with  gravimetric  moisture  content  0%,  5%, 
and  10%,  respectively. 


Table  Cl  Moisture  Contents  of  Soils  from  Site  1 


Depth  (inch) 

%  Water  by  Weight 

0 

3.1 

6 

K  4.6 

12 

4.6 

18 

4.3 

24 

4.4 

30 

4.6 

36 

4.9 
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It  can  be  seen  that  the  moisture  content  of  the  soil  is  approximately  5%,  although 
the  surface  soil  is  notably  drier.  Also,  it  was  noted  in  the  first  report  [1]  that  at  higher 
frequencies,  i.e.,  for  f  >  1  GHz,  the  higher  soil  absorption  decreases  the  penetration  depth, 
making  the  radar  insensitive  to  the  roots.  As  a  result,  the  data  from  Site  1  for  frequencies 
less  than  1  GHz  with  moisture  content  5%  are  of  special  interest  and  are  shown  in  Figures 
C5  and  C6.  From  the  last  two  figures,  the  average  values  of  dielectric  constant  and 
conductivity  at  six  inch  depths  are  found  to  be  approximately  5  and  0.0035  S/m, 
respectively.  These  soil  parameters  are  used  throughout  the  report. 
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Figure  C5.  Dielectric  constant  of  the  soil  samples  with  moisture  content  5%  from 
different  depths. 
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Figure  C6.  Conductivity  of  the  soil  samples  with  moisture  content  5%  from  different 
depths. 
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